{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Homogeneous Gas Reactor (Ignition Zero D problem)\n",
    "\n",
    "***Energy equation***\n",
    "$$\n",
    "\\frac{dT}{dt}= -\\frac{1}{\\rho c_p}\\sum_{k=1}^{N_{spec}}\\dot{\\omega_{k}} W_k h_k = S_T\n",
    "$$\n",
    "***Species equation***\n",
    "$$\n",
    "\\frac{dY_k}{dt}=\\frac{1}{\\rho}\\dot{\\omega_{k}}W_k=S_{Y_k},\\,\\,\\,k=1\\ldots N_{spec}\n",
    "$$\n",
    "\n",
    "## Environment Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "TChem_install_directory ='where/tchem/is/installed'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append(TChem_install_directory+'/lib')\n",
    "import pytchem\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import time\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.path.append(TChem_install_directory+'/example/runs/scripts/')\n",
    "import pmixSample\n",
    "import helper"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pre-simulation\n",
    "* Set variables; temperature and stoichiometric ratio (fuel/air).\n",
    "* Convert from stoichiometric ratio to mass fraction of CH4, O2, N2 and AR.\n",
    "* Create samples spaning over the variable ranges."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pressure, Temperature, phi(stoichiometric ratio)\n",
    "one_atm = 101325 # [Pa]\n",
    "Temp = 1300   #  temperature [K]\n",
    "Pressure = 1*one_atm; # [Pa]\n",
    "phi = 1.0; # Maximum equivalence ratio [-]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "N=100\n",
    "Nvar = 6\n",
    "sample = np.zeros([N,Nvar])\n",
    "fuel =\"CH4\"\n",
    "nC=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample[:,0] = Temp\n",
    "sample[:,1] = Pressure\n",
    "Yp_fuel, Yr_o2, Yr_n2, Yr_ar = pmixSample.getMassFraction(nC,phi)\n",
    "sample[:,2] = Yp_fuel\n",
    "sample[:,3] = Yr_o2\n",
    "sample[:,4] = Yr_n2 \n",
    "sample[:,5] = Yr_ar"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## TChem Simulation\n",
    "\n",
    "### Initialize TChem Driver Object\n",
    "\n",
    "* Initialization of Kokkos.\n",
    "* Create a TChem driver object. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "pytchem.initialize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "tchem = pytchem.TChemDriver()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Get help from TChem driver object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on TChemDriver in module pytchem object:\n",
      "\n",
      "class TChemDriver(pybind11_builtins.pybind11_object)\n",
      " |  A class to manage data movement between numpy to kokkos views in TChem::Driver object\n",
      " |  \n",
      " |  Method resolution order:\n",
      " |      TChemDriver\n",
      " |      pybind11_builtins.pybind11_object\n",
      " |      builtins.object\n",
      " |  \n",
      " |  Methods defined here:\n",
      " |  \n",
      " |  __init__(...)\n",
      " |      __init__(self: pytchem.TChemDriver) -> None\n",
      " |  \n",
      " |  cloneGasKineticModel(...)\n",
      " |      cloneGasKineticModel(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Internally create clones of the kinetic model\n",
      " |  \n",
      " |  computeGasEnthapyMass(...)\n",
      " |      computeGasEnthapyMass(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Compute enthalpy mass and mixture enthalpy\n",
      " |  \n",
      " |  computeGasNetProductionRatePerMass(...)\n",
      " |      computeGasNetProductionRatePerMass(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Compute net production rate\n",
      " |  \n",
      " |  computeGasReactionRateConstants(...)\n",
      " |      computeGasReactionRateConstants(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Compute forward/reverse rate constant\n",
      " |  \n",
      " |  computeJacobianHomogeneousGasReactor(...)\n",
      " |      computeJacobianHomogeneousGasReactor(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Compute Jacobian matrix for homogeneous gas reactor\n",
      " |  \n",
      " |  computeRHS_HomogeneousGasReactor(...)\n",
      " |      computeRHS_HomogeneousGasReactor(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Compute RHS for homogeneous gas reactor\n",
      " |  \n",
      " |  computeTimeAdvanceHomogeneousGasReactor(...)\n",
      " |      computeTimeAdvanceHomogeneousGasReactor(self: pytchem.TChemDriver) -> float\n",
      " |      \n",
      " |      Compute Time Advance for a Homogeneous-Gas Reactor\n",
      " |  \n",
      " |  createAllViews(...)\n",
      " |      createAllViews(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Allocate all necessary workspace for this driver\n",
      " |  \n",
      " |  createGasEnthapyMass(...)\n",
      " |      createGasEnthapyMass(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Allocate memory for enthalpy mass  (# samples, # species )\n",
      " |  \n",
      " |  createGasKineticModel(...)\n",
      " |      createGasKineticModel(self: pytchem.TChemDriver, chemkin_input: str, thermo_data: str) -> None\n",
      " |      \n",
      " |      Create a kinetic model from CHEMKIN input files\n",
      " |  \n",
      " |  createGasKineticModelConstData(...)\n",
      " |      createGasKineticModelConstData(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Internally construct const object of the kinetic model and load them to device\n",
      " |  \n",
      " |  createGasKineticModelConstDataWithArreniusForwardParameters(...)\n",
      " |      createGasKineticModelConstDataWithArreniusForwardParameters(self: pytchem.TChemDriver, reac_indices: numpy.ndarray[numpy.int32], factors: numpy.ndarray[numpy.float64]) -> None\n",
      " |      \n",
      " |      Creates clones of the kinetic model; modifies the Arrhenius forward parameters of the clones;and creates a const object of the kinetics models.factors is a 3d array of size: number of samples, number of reactions to be modified, and 3 kinetic parameters. kinetic parameters: pre exponential (0), temperature coefficient(1), and activation energy(2)\n",
      " |  \n",
      " |  createGasNetProductionRatePerMass(...)\n",
      " |      createGasNetProductionRatePerMass(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Allocate memory for net production rate per mass (# samples, # species)\n",
      " |  \n",
      " |  createGasReactionRateConstants(...)\n",
      " |      createGasReactionRateConstants(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Allocate memory for forward/reverse rate constants  (# samples, # reactions )\n",
      " |  \n",
      " |  createJacobianHomogeneousGasReactor(...)\n",
      " |      createJacobianHomogeneousGasReactor(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Allocate memory for homogeneous-gas-reactor Jacobian  (# samples, # species + 1  # species + 1)\n",
      " |  \n",
      " |  createRHS_HomogeneousGasReactor(...)\n",
      " |      createRHS_HomogeneousGasReactor(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Allocate memory for homogeneous-gas-reactor RHS  (# samples, # species + 1 )\n",
      " |  \n",
      " |  createStateVector(...)\n",
      " |      createStateVector(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Allocate memory for state vector (# samples, state vector length)\n",
      " |  \n",
      " |  freeAllViews(...)\n",
      " |      freeAllViews(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Free all necessary workspace for this driver\n",
      " |  \n",
      " |  getGasArrheniusForwardParameter(...)\n",
      " |      getGasArrheniusForwardParameter(*args, **kwargs)\n",
      " |      Overloaded function.\n",
      " |      \n",
      " |      1. getGasArrheniusForwardParameter(self: pytchem.TChemDriver, reac_indices: numpy.ndarray[numpy.int32], param_index: int) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrive pre exponential for reactions listed by reaction_indices\n",
      " |      \n",
      " |      2. getGasArrheniusForwardParameter(self: pytchem.TChemDriver, imodel: int, reac_indices: numpy.ndarray[numpy.int32], param_index: int) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrive pre exponential for reactions listed by reaction_indices\n",
      " |  \n",
      " |  getGasEnthapyMass(...)\n",
      " |      getGasEnthapyMass(self: pytchem.TChemDriver) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrive enthalpy mass per species for all samples\n",
      " |  \n",
      " |  getGasEnthapyMixMass(...)\n",
      " |      getGasEnthapyMixMass(self: pytchem.TChemDriver) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrieve mixture enthalpy for all samples\n",
      " |  \n",
      " |  getGasForwardReactionRateConstants(...)\n",
      " |      getGasForwardReactionRateConstants(*args, **kwargs)\n",
      " |      Overloaded function.\n",
      " |      \n",
      " |      1. getGasForwardReactionRateConstants(self: pytchem.TChemDriver, sample_index: int) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrive forward rate constants for a single sample\n",
      " |      \n",
      " |      2. getGasForwardReactionRateConstants(self: pytchem.TChemDriver) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrieve forward rate constants  for all samples\n",
      " |  \n",
      " |  getGasNetProductionRatePerMass(...)\n",
      " |      getGasNetProductionRatePerMass(*args, **kwargs)\n",
      " |      Overloaded function.\n",
      " |      \n",
      " |      1. getGasNetProductionRatePerMass(self: pytchem.TChemDriver, sample_index: int) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrive net production rate for a single sample\n",
      " |      \n",
      " |      2. getGasNetProductionRatePerMass(self: pytchem.TChemDriver) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrieve net production rate for all samples\n",
      " |  \n",
      " |  getGasReverseReactionRateConstants(...)\n",
      " |      getGasReverseReactionRateConstants(*args, **kwargs)\n",
      " |      Overloaded function.\n",
      " |      \n",
      " |      1. getGasReverseReactionRateConstants(self: pytchem.TChemDriver, sample_index: int) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrive reverse rate constants for a single sample\n",
      " |      \n",
      " |      2. getGasReverseReactionRateConstants(self: pytchem.TChemDriver) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrieve reverse rate constants  for all samples\n",
      " |  \n",
      " |  getJacobianHomogeneousGasReactor(...)\n",
      " |      getJacobianHomogeneousGasReactor(self: pytchem.TChemDriver, sample_index: int) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrive homogeneous-gas-reactor Jacobian for a single sample\n",
      " |  \n",
      " |  getLengthOfStateVector(...)\n",
      " |      getLengthOfStateVector(self: pytchem.TChemDriver) -> int\n",
      " |      \n",
      " |      Get the size of state vector i.e., rho, P, T, Y_{0-Nspec-1}\n",
      " |  \n",
      " |  getNumberOfReactions(...)\n",
      " |      getNumberOfReactions(self: pytchem.TChemDriver) -> int\n",
      " |      \n",
      " |      Get the number of reactions registered in the kinetic model\n",
      " |  \n",
      " |  getNumberOfSamples(...)\n",
      " |      getNumberOfSamples(self: pytchem.TChemDriver) -> int\n",
      " |      \n",
      " |      Get the number of samples which is currently used in the driver\n",
      " |  \n",
      " |  getNumberOfSpecies(...)\n",
      " |      getNumberOfSpecies(self: pytchem.TChemDriver) -> int\n",
      " |      \n",
      " |      Get the number of species registered in the kinetic model\n",
      " |  \n",
      " |  getRHS_HomogeneousGasReactor(...)\n",
      " |      getRHS_HomogeneousGasReactor(*args, **kwargs)\n",
      " |      Overloaded function.\n",
      " |      \n",
      " |      1. getRHS_HomogeneousGasReactor(self: pytchem.TChemDriver, sample_index: int) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrive homogeneous-gas-reactor RHS for a single sample\n",
      " |      \n",
      " |      2. getRHS_HomogeneousGasReactor(self: pytchem.TChemDriver) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrieve homogeneous-gas-reactor RHS_ for all samples\n",
      " |  \n",
      " |  getSpeciesIndex(...)\n",
      " |      getSpeciesIndex(self: pytchem.TChemDriver, species_name: str) -> int\n",
      " |      \n",
      " |      Get species index\n",
      " |  \n",
      " |  getStateVariableIndex(...)\n",
      " |      getStateVariableIndex(self: pytchem.TChemDriver, var_name: str) -> int\n",
      " |      \n",
      " |      Get state variable index\n",
      " |  \n",
      " |  getStateVector(...)\n",
      " |      getStateVector(*args, **kwargs)\n",
      " |      Overloaded function.\n",
      " |      \n",
      " |      1. getStateVector(self: pytchem.TChemDriver, sample_index: int) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrieve state vector for a single sample\n",
      " |      \n",
      " |      2. getStateVector(self: pytchem.TChemDriver) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrieve state vector for all samples\n",
      " |  \n",
      " |  getTimeStep(...)\n",
      " |      getTimeStep(self: pytchem.TChemDriver) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrieve time line of all samples\n",
      " |  \n",
      " |  getTimeStepSize(...)\n",
      " |      getTimeStepSize(self: pytchem.TChemDriver) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrieve time step sizes of all samples\n",
      " |  \n",
      " |  modifyGasArrheniusForwardParameters(...)\n",
      " |      modifyGasArrheniusForwardParameters(self: pytchem.TChemDriver, reac_indices: numpy.ndarray[numpy.int32], factors: numpy.ndarray[numpy.float64]) -> None\n",
      " |      \n",
      " |      Modify the cloned kinetic models Arrhenius parameters\n",
      " |  \n",
      " |  setNumberOfSamples(...)\n",
      " |      setNumberOfSamples(self: pytchem.TChemDriver, number_of_samples: int) -> None\n",
      " |      \n",
      " |      Set the number of samples; this is used for Kokkos view allocation\n",
      " |  \n",
      " |  setStateVector(...)\n",
      " |      setStateVector(*args, **kwargs)\n",
      " |      Overloaded function.\n",
      " |      \n",
      " |      1. setStateVector(self: pytchem.TChemDriver, sample_index: int, 1d_state_vector: numpy.ndarray[numpy.float64]) -> None\n",
      " |      \n",
      " |      Overwrite state vector for a single sample\n",
      " |      \n",
      " |      2. setStateVector(self: pytchem.TChemDriver, 2d_state_vector: numpy.ndarray[numpy.float64]) -> None\n",
      " |      \n",
      " |      Overwrite state vector for all samples\n",
      " |  \n",
      " |  setTimeAdvanceHomogeneousGasReactor(...)\n",
      " |      setTimeAdvanceHomogeneousGasReactor(self: pytchem.TChemDriver, tbeg: float, tend: float, dtmin: float, dtmax: float, jacobian_interval: int, max_num_newton_iterations: int, num_time_iterations_per_interval: int, atol_newton: float, rtol_newton: float, atol_time: float, rtol_time: float) -> None\n",
      " |      \n",
      " |      Set time advance object for homogeneous gas reactor\n",
      " |  \n",
      " |  showViewStatus(...)\n",
      " |      showViewStatus(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Print member variable view status\n",
      " |  \n",
      " |  ----------------------------------------------------------------------\n",
      " |  Static methods inherited from pybind11_builtins.pybind11_object:\n",
      " |  \n",
      " |  __new__(*args, **kwargs) from pybind11_builtins.pybind11_type\n",
      " |      Create and return a new object.  See help(type) for accurate signature.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(tchem)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create Kinetic Model \n",
    "\n",
    "* Inputs are the reactions mechanism files; in this case, we use the GRI3.0 gas reaction mechanism"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "inputs_directory = TChem_install_directory + '/example/data/ignition-zero-d/gri3.0/'\n",
    "tchem.createGasKineticModel(inputs_directory+'chem.inp',inputs_directory+'therm.dat')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "tchem.setNumberOfSamples(N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ### Model Variation\n",
    "* Set list of reaction to be modified ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nparameter= 10\n",
    "reaction_index = np.zeros(Nparameter,dtype=int)\n",
    "for i in range(Nparameter):\n",
    "    reaction_index[i]=i + 3  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Set a 3d array, factors, of dimension: Number of samples (N), number of reaction to be modified (Nparameter), 3 that corresponds to three Arrenius parameters: pre exponetial (0), temperature coeficient(1), activation energy(2). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "μ=0 \n",
    "σ=.25\n",
    "A_modifier = np.random.lognormal(μ, σ, size=(N, Nparameter))\n",
    "factors = np.empty([N,Nparameter,3])\n",
    "ones = np.ones([N,Nparameter])\n",
    "factors[:,:,0] = A_modifier # //pre exponetial \n",
    "factors[:,:,1] = ones # temperature coef \n",
    "factors[:,:,2] = ones # activation energy "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Create a list of const objects of the kinetic model with the modified parameters. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "tchem.createGasKineticModelConstDataWithArreniusForwardParameters(reaction_index,factors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Verify the changes comparing reference values vs modified variables for a specific sample (isample)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "ver_reaction_index = np.zeros(3,dtype=int)\n",
    "ver_reaction_index[0] = 1\n",
    "ver_reaction_index[1] = 3\n",
    "ver_reaction_index[2] = 6\n",
    "pre_exponetial_index = 0\n",
    "isample=84\n",
    "\n",
    "A_original = tchem.getGasArrheniusForwardParameter(ver_reaction_index, pre_exponetial_index); \n",
    "A_modified = tchem.getGasArrheniusForwardParameter(isample, ver_reaction_index, pre_exponetial_index);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pre exponential Original: [5.e+17 2.e+13 8.e+13]\n",
      "Pre exponential Modified: [5.0000000e+17 2.5474552e+13 6.5018563e+13]\n",
      "Factor                  : [1.         1.2737276  0.81273204]\n"
     ]
    }
   ],
   "source": [
    "print('Pre exponential Original:', A_original)\n",
    "print('Pre exponential Modified:', A_modified)\n",
    "print('Factor                  :', A_modified/A_original)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "for isample in [3, 20, 40, 53, 83]:\n",
    "    A_original = tchem.getGasArrheniusForwardParameter(reaction_index, pre_exponetial_index); \n",
    "    A_modified = tchem.getGasArrheniusForwardParameter(isample, reaction_index, pre_exponetial_index);\n",
    "    factor = A_modified/A_original\n",
    "    for i, ifac in enumerate(factor):\n",
    "        error = (ifac - A_modifier[isample,i])/ A_modifier[isample,i]\n",
    "        if abs(error) > 1e-15:\n",
    "            print('error in ',error,' in sample No',isample,' factor No', i)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set State Vector\n",
    "\n",
    "* Get index for variables. \n",
    "* Pass a numpy array to the TChem object to set the state vector. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "Variables = ['T','P','CH4','O2','N2','AR']\n",
    "indx=[]\n",
    "for var in Variables:\n",
    "    indx += [tchem.getStateVariableIndex(var)]\n",
    "\n",
    "state = np.zeros([N, tchem.getLengthOfStateVector()])\n",
    "for sp in range(N):\n",
    "    state[sp,indx] = sample[sp,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "tchem.createStateVector()\n",
    "tchem.setStateVector(state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set Time Integrator and Its Parameters \n",
    "* Kokkos team policy is constructed for parallel execution\n",
    "* Workspace per team is also allocated with the teampolicy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "tend = 0.05\n",
    "tchem.setTimeAdvanceHomogeneousGasReactor(tbeg=0, \n",
    "                     tend=tend, \n",
    "                     dtmin=1e-10,\n",
    "                     dtmax=1e-3, \n",
    "                     atol_time=1e-12,\n",
    "                     rtol_time=1e-6,\n",
    "                     max_num_newton_iterations=20,\n",
    "                     atol_newton=1e-18,\n",
    "                     rtol_newton=1e-8,\n",
    "                     num_time_iterations_per_interval=20,\n",
    "                     jacobian_interval=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Advance simulations: for each iteration, get state vector. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time is 1.6960 mins for time integration\n"
     ]
    }
   ],
   "source": [
    "solution = []\n",
    "time_simulation = []\n",
    "\n",
    "tavg = 0\n",
    "t0 = time.time()\n",
    "icount = 0\n",
    "while (tavg < tend*0.999):\n",
    "    tavg = tchem.computeTimeAdvanceHomogeneousGasReactor()\n",
    "    solution += [tchem.getStateVector()]\n",
    "    time_simulation += [tchem.getTimeStep()]\n",
    "    icount +=1\n",
    "    if (icount%50 == 0):\n",
    "        print(f\"count {icount:3d} simulation time {tavg:10.3e} s, total time {tend:10.3e} s\")\n",
    "\n",
    "t1 = time.time()\n",
    "\n",
    "print(f\"Wall time is {(t1 - t0)/60:0.4f} mins for time integration\")\n",
    "\n",
    "solution = np.array(solution)\n",
    "time_simulation = np.array(time_simulation) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Post-processing\n",
    "* Get variable index for plotting. \n",
    "* Compute the ignition delay time with data produced by the simulation.\n",
    "* Plot time profiles for important variables and for the ignition delay time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "T_indx   = tchem.getStateVariableIndex('T')\n",
    "CH4_indx = tchem.getStateVariableIndex('CH4')\n",
    "O2_indx  = tchem.getStateVariableIndex('O2')\n",
    "CO_indx  = tchem.getStateVariableIndex('CO')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "T_threshold=1500\n",
    "def computeIgnitionDelayTime(T_threshold=1500):\n",
    "    ntimes, nsamples, nVars = np.shape(solution)\n",
    "    ignDelayTime = []\n",
    "    for isample in range(nsamples):\n",
    "        Indx = np.where(solution[:,isample,T_indx] >= T_threshold )\n",
    "        ignDelayTime += [time_simulation[Indx[0][0],isample] ] \n",
    "    return np.array(ignDelayTime)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "IgnTime = computeIgnitionDelayTime()\n",
    "np.savetxt('QoI_GRI3.txt',IgnTime)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmsAAAEcCAYAAACYg/MAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABjmElEQVR4nO3dd3hUZfbA8e/JpAAJ0pGqoamoICKiQYSIIoK61p9lcRXrIjYsu3bFsmJl0bVXXNfelaIiEFAMKiIqiApKkN5rejLn98c7k0wmk2SSTDIp5/M895mZe9+5895Lwpy85byiqhhjjDHGmLopJtoVMMYYY4wxZbNgzRhjjDGmDrNgzRhjjDGmDrNgzRhjjDGmDrNgzRhjjDGmDouNdgXqopiYGG3atGm0q2GMMQ1SVlaWqqo1FhgTJgvWQmjatCmZmZnRroYxxjRIIpId7ToYU5/YXzbGGGOMMXWYBWvGGGOMMXWYBWvGGGOMMXWYBWvGGGOMMXWYBWvGGGOMMXWYBWvGGGOMMXWYpe4wUZGZmck333yDqhIXF1dii42NLbUveL/H44n2JRhjjDG1woI1U2uysrKYPn06b731FlOnTiU7u+qplkSk0gFeVYLC6hwTkUpdTyTL1VTZ+nLO+vT5xhhTEVHVaNehzonpEKNdb+rqnksMgiAixEgMHvHgifGEfIyNiSXOE0e8J564GN9jwOtQ+/yvm8Y2pVlcMxLjE91jXGKZr2Nj6k+MnZ2dzSeffMJbb73Fxx9/TGZmJu3bt+eMM87g5JNPplmzZuTn54fcCgoKanR/OO8pLCyM9i00piHKUtXEaFfCmPqi/nzr16KYghiGdRuGqqIoXvWi6h4LtZBCb2HIxwJvAbkFuezJ20NeYR75hfnu0Zsf8nWBt6BK9Yv3xBcFcEnxSbRPbE/7xPbsnbi3e0zau9Tr5vHNa/0v/XXr1nHwwQezfft22rZty3nnncdZZ53FkCFDiI2tHz96Xq+3REBXmYAwXOH+wVSZP6xqomx9OWd9+vzG6vTTT492FYypV6xlLYTExEStjeWmVJUCbwF5hXlkF2STlZ9FZl6me8zPDOv1rtxdbM7azMY9G9mUuYmt2VtDflaT2CYlA7jEvdmvzX4M6DSAwzodRssmLSN+fTNmzGDUqFE899xzjBkzpt4EaMaYmiUi1rJmTCXYt2cUiQhxnjjiPHEkxkfm/638wny2ZG1hY6YL3vxB3MbMjUX71u1ex6L1i3hx8YtF7+vVuhcDOg0o2vp37E9SfFK16rJq1SoARo4caYGaMcYYU0X2DdrAxHni6Ni8Ix2bd6yw7LbsbXy37jsWrlvIwvUL+fLPL3l9yesACELvdr0Z3HUw16ZcywFtD6h0XTIyMoiLi6Njx4rrYowxxpjQrBs0hNrqBq2LNu7ZyHfrfQHcuoXMWjmL7Pxszjn4HG4fcju92/UO+1znnnsu3377LStWrKjBGhtj6hvrBjWmcixYC6ExB2vBNmVu4pGvHuGJb58gKz+Lsw46i9uH3M5B7Q+q8L0pKSkkJiby+eef10JNjTH1hQVrxlSOrWBgytU+sT0PDH+AjPEZ3DT4JqYtn0afp/pwzjvnsDt3d7nvXbVqFfvuu28t1dQYY4xpmCxYM2Fp26wt9x17HxnXZHDL0bfwzs/vMObDMWWmKMjJyWH9+vUWrBljjDHVZMGaqZQ2zdpw77B7eXD4g7y37D0enP9gyHKrV68GIDk5uRZrZ4wxxjQ8FqyZKrn2yGs5+6CzuWX2Lcz8fWap4/60HdayZowxxlSPBWumSkSE5//yPL3b9uacd89h1Y5VJY5nZGQA1rJmjDHGVJcFa6bKkuKTeP/s98nMy2Tygskljq1atQqPx0Pnzp2jUzljjDGmgbBgzVRLrza9OLb7sXz464clJhtkZGTQpUsXW7nAGGOMqSYL1ky1nbr/qazcsZIlm5YU7bO0HcYYY0xkWLBmqu3k/U9GED745YOifRkZGRasGWOqRUTGichKEckRke9E5OhyyqaKyIcisl5EskTkRxG5KES5ob5z5YjIHyIytmavwpjqs2DNVFuHpA4c0eUIPvz1QwDy8/NZu3atTS4wxlSZiJwNPArcBxwKfAXMEJF9ynjLIOAn4EzgYOAp4FkR+WvAObsB033nOhSYCPxHRM6oqeswJhIsWDMRcer+p/Ld+u9Yu2sta9aswev1WsuaMaY6rgOmqOpzqrpMVa8C1gOXhyqsqvep6m2qOl9V/1DVp4D3gMBAbCywTlWv8p3zOeBl4IYavhZjqqVWgzURuVlEvhWRXSKyWUQ+FpGDg8pMEREN2hYElUkQkf+IyBYRyRSRj0SkS1CZViLyiojs9G2viEjLWrjMRunILkcCsGzLMv7880/AcqwZY6pGROKBw4DPgg59hmtBC9dewPaA1ykhzvkpMEBE4ipbT2NqS21P1UsFngS+BQS4G/hcRA5U1W0B5T4H/hbwOi/oPJOBU4Bzga3AJGCqiBymqoW+Mq8B+wAjAQWeB14BTo7g9Rif5JbJAGTsyGDvXXsD0KpVqyjWyBhTh8WKyMKA18+q6rMBr9sCHmBj0Ps2AseF8wEichJwLHBUwO4OuO+X4HPG+j5zfTjnjoSYmBht2rRpbX2cqaOysrJUVStsOKvVYE1VRwS+FpG/ATtxv0wfBxzKVdUNoc4hIi2Ai4ELVXVmwHlW4X6JPxWR3sAJwGBV/cpX5u/AFyKyv6r+GtkrM5336oxHPGTsyKBFTgsAmjRpEuVaGWPqqAJVHRBGueDFhyXEvlJE5CjcH+xXq+o3YZwz1P4a1bRpUzIzM2vzI00dJCLZ4ZSL9pi15r46bA/aP1hENonIbyLynIi0Dzh2GBBHQFO2qq4GllHcPJ4C7MENIvWbD2RSuSZ0E6bYmFi6tujKyh0rycnJASxYM8ZU2RagENcSFqg9pVvbShCRwcAM4A7fuLVAG8o4ZwGul8aYOinawdqjwGIgPWDfJ8D5uObr64GBwGwRSfAd74D7Jd4SdK6NFP8SdgA2a0CWVt/zTZT+RQVARC4TkYUisrCgoKA619RoJbdMJmNHhgVrxphqUdU84DtgeNCh4ZT8I7wEERmCC9TuUtXJIYqkU7obdTiwUFXzq1xhY2pY1II1EZkEDAbOCBhnhqq+oaofqepPqvoxbszZ/sCJFZ2Sks3YoZq0y2xCV9VnVXWAqg6wrPtV061lNwvWjDGRMgkYIyKXiEhvEXkU6AQ8DSAiE0Vklr+wiKTiArWngVdFpINvaxdwzqeBLiIy2XfOS4AxwMO1ckXGVFFUohIR+TdwDnCMqv5RXllVXScia4Bevl0bcANP2wKbA4q2B+YFlGkvIuJvXRMRAdpRQRO6qbrklsms272O3U13AxasGWOqTlXfFJE2wG1AR2AJMEpVV/mKdAR6BLxlDNAMl4YjMBXHKiDZd86VIjIK+DcuBcg63Li2d2vuSoyB9PR00tLSGDp0KEcccQSFhYV4vd6w31/rwZrvr6NzgFRV/SWM8m2BzhTP0vkOyMc1Xb/mK9MF6E1x83g6kIQbu+bflwIkUk4TuqmejkkdAdiW4yb2JiQklFfcGGPKpapP4jIIhDo2JsTrMaHKBpWbC/Svfu2MKV9BQQHp6ek899xz/O9//yuxfnZl1WqwJiJP4FJynApsFxH/+LE9qrpHRJKACcC7uOAsGZdhehPwPoCq7hSRF4CHRGQTxak7fsQ3JVtVl4nIJ8AzInIprvvzGWCqzQStOXsl7AXArtxdxMbG2iLuxhhjGpXNmzfzySefMG3aND799FN27NhBTExMUaAmIgwbNoxjjjkGj8fDzTffHNZ5a/vbdJzvcVbQ/rtwQVoh0Ac3waAlLmCbA5ylqrsDyl+Lm73zJtDUd77zA8e+AaOBxyieNfoRcGWErsOE0DyhOQC783dbF6gxxpgGz+v1smjRIqZPn860adP49ttvUVX23ntvTjvtNEaNGkWLFi045ZRTyMvLIz4+nnvuuYeUlBSAuhmsqapUcDwbGFFeGV+5HOAq31ZWmW3AeZWto6k6f8taZl6mBWvGGFMOb4yX1TtX0yyuGc3imtEktgluaHXZ0lenk5aRRmpyKildU2qppibYzp07+eyzz5g+fTozZsxg48aNiAgDBw7krrvuYtSoURx66KHExBTP4Zw1axZpaWmkpqYWBWqVYf1UJmKKgrUCC9aMMaY8Oc1z2GdyyTXp/YFbqC2nIIf5f87Hq148MR6uS7mOwV0H06l5Jzo270j7xPbExthXek1QVX7++eei1rP58+dTUFBAy5YtOeGEExg1ahQnnHAC7dq1K/McKSkpVQrS/Oxf1kRMUbBWaMGaMcaUJz4rnidPfpKs/KzSW0HJ17tyd5GxI4NC30ifAm8BD85/kAd5sOh8MRJD+8T2dGreyQVwSR1LPjZ3jxbUhScrK4vZs2czffp0pk+fzqpVbhJy3759+cc//sGoUaM48sgja21stv2LmYjxB2vZhdkWrBljTDli82O5uP/FYZdPX53Osf89lrzCPOI98bx+xut0bN6R9bvXs273OtbvKX5cs2sN3679lk2Zm9Cg1KKCsHfS3mUGc/7Xeyftzbdrv21U3a5//PFHUevZnDlzyM3NJTExkeOOO45bbrmFUaNG0aVLl6jUzYI1EzHN490Eg2zNpkWTFlGujTHGNBwpXVOYdf6sSgVP+YX5bMrcVCKYW7d7nQvw9rjnC9ctDBnUBYqNieX+Y+/n0sMuLfqjvCHIy8vjiy++KArQfv3VJYvo1asXY8eO5cQTT2TIkCF1Ig2VVCfvR0OVmJiotsBu1TT9V1ParWxH8m/JzJs3r+I3GGMaHRHJUtXEaNcjmurS90yBt4CNezYWt87tXs9bS99iTsacEkFcjMRwaIdDGbrvUIYmD+XofY6mVdNWUax55a1bt44ZM2Ywbdo0Zs6cyZ49e4iPjyc1NZVRo0YxatQoevXqVfGJIiTc3wVrWTMRlRSfRJ7mWTeoMcbUE7ExsXTeqzOd9+pctK/v3n1LdLs+MPwBtmRuYe6quTzx7RNMWjAJQei7d9+i4G3IvkNo26xtFK+ktPT0dF555RWysrL44YcfWLx4MQBdunThr3/9KyeeeCLDhg0jKSkpuhWtgAVrJqISPAlkq41ZM8aY+qy8btecghy+WfsNczPmMnfVXJ5b9ByPffMYAAe1O4gh+w6hY1JHCrSAE3qcELXxbunp6QwdOpT8/HwADjnkECZOnMioUaPo06dPhalS6hLrBg2hLjVP1zfdH+3O1sVbGZE5grfeeiva1THG1EHWDdqwvmfyCvNYuG4hczPmMu/PeczNmEt2QTYAcTFxzL5gNoP3GVzr9TrrrLN4++23AfB4PNxzzz1hJ6GtLeH+LsRUVMCYykiITaCQQmtZM8aYRiLeE8+groO4+eibmTF6BrcefSsxvvAi35vP/739f6RlpNVqnWbNmsV7771HTEwMHo+naFxafWXBmomoeE88BVpgwZoxxjRSw7oNIyE2AY94iPfEg8IxLx/D6W+ezh/b/6jxz//xxx85/fTT6d27NzNmzOCee+5h1qxZ1UpKG202Zs1EVIIngUKxljVjjGmsgse79evQj0npk5j45UR6P9Gb8UeMZ3iP4Xy79tuI53BbvXo1o0aNonnz5kyfPp2uXbty/PHHR+z80WLBmomoeE88XrwWrBljTCOW0jWlRBB265BbufDQC7ll1i08+NWDPPjVg8QQQ0JsArPOnxWRgG3Hjh2MHDmS3bt38+WXX9K1a9dqn7OusG5QE1EJngS8MRasGWOMKalT805MOXUK4w4fB4AXL3mFeREZz5abm8tpp53Gb7/9xvvvv0+fPn2qfc66xII1E1GxMbHgwYI1Y4wxIZ3X5zyaxLrvCK966bt332qdz+v1cuGFF5KWlsZLL73EsGHDIlHNOsWCNRNRcRJnwZoxxpgypXRNYfb5s7ny8CuJ88Rx//z7ySvMq/L5brnlFl5//XUmTpzI6NGjI1jTuqPMMWsi8mcVzqfAiaq6pOpVMvWZB48Fa8YYY8rlH9M2eJ/BnPPuOZz19lkc0fmISk84eOKJJ3jggQe4/PLLufHGG2uwxtFV3gSDLsB0YHOY54oBzgPiq1spU48pEEOdWPjWGGNM3Xb2wWcz7bdpvPLTK3z060c0iW0S9oSDDz/8kKuvvpqTTz6Zxx57rF6tSFBZFc0GvVtVvwnnRCISC/yt+lUy9VmMxkAMxMXFRbsqxhhj6oH92+4PgKJFEw4qCtYWLFjAueeey4ABA3j99deJjW3YyS3KG7N2O7A63BOpaoHvPWurWylTfwkCQoP/xTHGGBMZw7oNI8HjemNiJIbU5NRyyy9fvpyTTz6ZTp068fHHH5OY2PBXLiszWFPVf6nq+nBOIiJJAe/ZGKnKmfonBteyZsGaMcaYcPgnHHRM6kjXFl05ssuRZZbdtGkTI0eOBOCTTz6hffv2tVXNqCozWBORi8M5gYi0BGZGqkKmnvONWbNgzRhjTLgG7TOIu4+5mz+2/8FXq78KWWb27Nn079+f1atXM3XqVHr27FnLtYye8rpBnxWRcufAikhrYBZQvSQppsEQtW5QY4wxlXfuweeSGJfIFdOvIH11eolj6enpHH/88axd60Zaeb3eaFQxasr7Rv0f8JKI5Kjqu8EHRaQdLlDrDpxSQ/Uz9Ywg1rJmjGlwRCQGGAjsA5TKTaSq/631SjUwP278kdzCXH7Y+APD/juM2efPLppo8NZbb1FYWAhAYWEhaWlp9Xph9soq7xt1DC4Nx2sicqaqfuw/ICIdcYFaF2CUqs6r0VqaesPfsubxeKJdFWOMiQgRORD4AOgBhMoPoYAFa9WUlpGGV12LWfCs0JUrVwLuuyU+Pp7U1NRoVTMqygzWVFVF5DzgLeAtETlFVT8TkX1wgVo7YISqppd1DtP4iLqWNQvWjDENyJO478uzgJ+A3OhWp2FKTU4lwZNAdkF2iVmh27Zt47PPPuOkk05i0KBBpKamNqpWNaggz5qqForI2cD7wPsicgUwAWgOHKuq39V8FU194g/WrBvUGNOA9AfGqOp70a5IQ5bSNYVZ58/ijLfOoPNenYta1W677Tays7M5++yzOe+886Jcy+iocG1QX/60M4AvgBeApkCqBWqmTNYNaoxpWLYAVV+80oQtpWsKf9n/L/y29Te86mX+/Pk8/fTTAFx22WWkpzfOzrzy1gYN7n/fjeur/wP4R9CyDqqqF0S+eqa+sTFrxpgG6N/AFSIyQ1ULo12Zhu6IzkfwzHfP8OuWX7n11ltRVQDy8vIa3cQCv/L6qobgBk0GWgV08G2BgsuZxsqXZ82CNWNMA9IO2B/4WURmAtuCjquq3ln71WqYBnYeCMB5L5zHot8XFe2PiYlpdBML/MqbYJBci/UwDYUvbG/IC+oaYxqd2wKe9wpxXAEL1iJke852ABblLIILgJeBNdC8efNG2aoGFS/kbkzl+IK1GE+FwyGNMaZeUFX7D60WfbHqC/dEcCPrk4E1kJmZGb1KRZn9AJqIEl8KIrWecWOMMVWQmpxKjMS4P/69QIbbv88++0SxVtFV3tqghSIyMNwTiYjH957+kamaqZf8MZr1ghpjGhgROUlEHhKRF0TkQRE5Mdp1aohSuqYwsudIYr2xRV2gAMcee2xU6xVN5bWsCdBBRPYJZwP2xb6iGz2xH4FGTUTK3caMGRPtKkbclClTSEpKinY1TA0SkeYiMhf4CLgGGAWMBz4SkTQRsR+ACOvfsT+FnkJY5/5fiY+P5/zzz492taKmojFr71fyfNb31cj5Jxao2I9CY7R+/fqi51OnTuXSSy8tsa9p06bRqFaV5OXlER8f3+A/sy4TkXHAP4COwFJgvKp+UUbZJsDTuAS2vYH5qpoaVCYVmBPi7b1V9ZdyqnKf77x/A97wJYz3AOcAT/mOXx32hZkKdW7e2Q2nSYJDuh/C5Zdf3mgnF0D5LWsXAhdVYVtZg/U1dZy/Zc2/vptpXDp06FC0tWzZstS+efPmcdhhh9GkSRO6devGrbfeSl5eca7R5ORk7r77bsaMGUPz5s3p2rUrb775Jjt27OCcc84hKSmJXr168dlnnxW9Jy0tDRFh6tSp9OvXjyZNmnDYYYfx3Xcl83Z/9dVXDB06lGbNmtG5c2cuv/xydu3aVXQ8NTWVyy+/nBtuuIF27dpx1FFHATBp0iT69u1LYmIinTt35pJLLmHHjh1Fn33hhReSmZlZ1Ho4YcKEomt5+OGHS9QhNTWVK6+8ssT1TpgwgYsuuoiWLVsyevTosOraGPhWz3kUFwgdCnwFzPD15ITiAXKAx4FpFZz+IFwA6N+WV1D+DOA2VX3Vn2dNVQtV9VXgdt9xE0Fd9urinhwNP2z7gfHjxzfahLhQTrCmqi9Xcdtemxdg6paiCQZqLWumpE8//ZTRo0dz5ZVXsnTpUl588UXeeecdbrnllhLlJk+ezMCBA1m0aBFnnXUWF1xwAX/9618ZNWoUixcvZsiQIZx33nnk5OSUeN8NN9zAAw88wMKFC+nevTsnnngiWVlZAPz0008cf/zx/OUvf+GHH37gvffeY/HixVx00UUlzvG///0PVeWLL77gv/91ecFjYmKYPHkyS5cu5bXXXuObb77hqquuAmDQoEFMnjyZZs2asX79etavX88NN9xQqfsyadIkDjjgABYuXMh9990Xdl0bgeuAKar6nKouU9WrgPXA5aEKq2qmqo5V1WcpGuVUpk2quiFgqyjRbRvg5zKO/ew7biJoS9YW9+Qw0L8pue1ySUtLi2qdoslSd5iIstmgNWv8+PEsXry4Vj+zX79+TJ48udrn+de//sU//vEPLrzwQgB69OjBAw88wHnnncdDDz1U1IU+YsQIxo0bB8Bdd93FpEmT6NmzZ9F4ldtvv50XX3yRJUuWMGDAgKLz33777YwYMQKAl156iS5duvDaa69xySWX8NBDD3H22Wdz/fXXF5V/6qmnOPTQQ9m0aRPt27cHoFu3bjzyyCMl6j1+/Pii58nJyTz44IOccsopvPzyy8THx9OiRQtEhA4dgnOFh2fo0KH885//LHp9/vnnh1XXhkxE4oHDgIeDDn0GDIrARywUkQRcoHWvqobqGg20EjgJmBni2CisRynift36q3sSAyjE9Gi8CXHBgjUTYTHiGmutG9QE++677/jmm2944IEHivZ5vV6ys7PZsGEDHTt2BKBv375Fx5OSkmjWrBl9+vQp2rf33nsDsGnTphLnDxzPkpSURJ8+ffj555+LPnvFihW8+eabRWX8rb+///57UQB02GGHlar37NmzmThxIsuWLWPnzp0UFhaSl5fHhg0b6NSpU9VuRoDAgLMyda3nYkVkYcDrZ30tYn5tcd2aG4PetxE4rhqf62+Z+xaIx41BmyUiqao6r5z3PQM84ptI8KrvPB1wY9YuwbUCmgg6vvvxTPxyInghNiaWJ/7xRKMes2bBmoks9T9Yy1pNiEQLV7R4vV7uvPNO/u///q/UsXbt2hU9j4uLK3FMRErs87fAeb3h/0Hg9Xq55JJLuPbaa0sd69y5c9HzxMTEEsdWrVrFiSeeyKWXXsrdd99NmzZtWLRoEeeee26JsXahxMTElBoOkJ+fX6pc8GeGW9d6rkBVB1RcrNR/JBJiX9hU9Vfg14Bd6SKSDNwAlBmsqeq/RaQdcC0wJqAuucD9qvpoVetkQhuaPBRB0FXKlIumMHrI6GhXKaosWDORZZk7TBn69+/PL7/8Qs+ePWvk/AsWLKB79+6Ay3S+ZMmSoq7T/v37s3Tp0kp/9sKFC8nLy+Pf//530Xq3U6dOLVEmPj6ewsLSQ57atWtXYiZsTk4Ov/zyC4ceemi5n1nVujYwW4BCSq9D3Z7SrW3V9TWuhaxcqnqLiDwEHAm0xq0PusDGadcMEaFZTDMyN2UysGPYKV8bLFvBwESUjVkzZbnjjjt47bXXuOOOO1iyZAm//PIL77zzTonxWtVx7733MnPmTJYuXcpFF11EfHw8f/3rXwG48cYb+eabbxg7dizff/89K1asYOrUqfz9738v95y9evXC6/UyefJkVq5cyeuvv16qdTM5OZmcnBxmzpzJli1biiY1DBs2jFdffZW0tLSiOoVqWQtW1bo2JKqaB3wHDA86NBw3KzSS+uG6NSukqttVdYZvVugMC9RqVmxhLCS4CUKNXaWCNRGJEZGDRWSoiCRW/A7T2FiwZsoyYsQIpk2bxpw5cxg4cCADBw7k/vvvj9gSMvfffz/XX389/fv3Z/ny5UydOrWoi7Fv377MmzePjIwMhg4dyiGHHMLNN99cNP6tLH379uXRRx9l0qRJHHjggTz//POl0nEMGjSIsWPHcu6559KuXTsefPBBAG6++WaGDRvGKaecwvHHH8/gwYPp37/iBV6qWtcGaBIwRkQuEZHeIvIo0AmXSw0RmSgiswLfICIHikg/3Ji3JBHp53vtPz5eRE4VkV4icpCITAROxaX7IOhcQ/zJbn3Py91q5hY0Xunp6ezctBMSYPTo0Y06bQeAhJtiQUSuAO6keIry4aq6SEQ+AGar6mM1U8Xal5iYqI15wdjquOTZS3hh/QvMGTWH1MNTo10d0wikpaVxzDHHsHnzZtq2bRvt6pgwiEiWqlb4B78vKe4/cbnQlgDX+icCiMgUIFVVkwPKZ+BW0ylBVcV3/J/AZUBnIBuXaHeiqk4P8dle4EhV/cb3vKwvS3EfoZ6KrieQfc+Ub+LEidyy4hYoAM+rHu655x5uvvnmaFcr4sL9XQhrzJqIXIpLTvgibur0WwGHv8AlBGwwwZqpOltuyhgTKar6JPBkGcfGhNiXXMH5HgQeDPPjj6E4t9owbIWeWpWamgp/AG3Ak+xp1Gk7IPwJBtcBj6jqjb4lNgL9glsOxJhiFrMZY+oxVZ0b8DwtilVpnLqAdBYURS4Q6BLtCkVXuGPWugGflnEsE2gZkdqYes9a1kxtS01NRVWtC9TUGBH5Q0QOKePYwSLyR23XqaFLy0hzY58FCrSAtIy0aFcpqsIN1rYAyWUc2x9YG5HamAbDlpsyxjQgyUBCGceaEGKcnKme1ORU98e/QrwnntTk1GhXKarCDdY+Bu4Qke4B+1RE2uKSBH4Q6YqZ+s1mgxpjGpiy/lMbAOyoxXo0CildU+i4vSMUwOT+k0np2nhXL4Dwg7XbcJmalwCf435oHwOW4RIX3l0jtTP1TlE3qPWGGmPqMRG5VkT+FJE/cd95H/tfB2ybgSeAT6Jb24YnPT2d9b+49Hfjzxjf6FN3hDXBQFW3isgAYDwwAvjd997HgX+r6q4aq6Gpl6wb1BhTz/0B+PO4XQAsBDYHlcnFzRh9vhbr1SikpaWheQpxkJuXS1pamq0NWh4RiQceAF5T1XuAe2q8VqbesqS4xpiGQFU/BD6EovVo71bVlVGtVCOSmpqKzHCzQeObxjf61B0VdoP6lv34O9C05qtjjDHG1Dl/BzaFOiAiiSISV8v1afBSUlJo2aUlAPe/dn+jblWD8MesfQ/0qcmKmAbGxqwZYxqO53xbKM/4NhNB6avT2d7TLb160483kb66cY9ZCzdYux64QUROEl97sDGhFE0wsF7QRm/jxo1cc8019OjRg4SEBDp37szIkSOZPt2t7JOcnFxqnU2Ahx9+mOTk5JDn/PLLL4mNjeXggw+uyaobE+wYfF2iIXwEHFuLdWkU0jLSiiKU/ML8Rp9nLdwVDN4GWuB+WAtEZBMlv45VVS3PjCn6qbAxa41bRkYGRx11FM2bN2fixIkccsgheL1eZs2axdixY/nzzz8rfc7t27dz/vnnc+yxx7J2raV2NLWqPWV0g+ImHexdi3VpFFKTU8ELeCDOE9fo86yFG6zNIgJtJSJyM3A6LpFuLrAAuFlVlwSUEdyC8ZcBrYCvgStUdWlAmQTgYeBc3Fi6WcA4VV0TUKYVLr3IX3y7PgKuUtUd1b0OUzFrgG3cxo0bh6qycOFCkpKSivb37t2b0aNHV+mcF198MRdccAGqyjvvvBOpqhoTjk24oUBzQhzrA2yt3eo0fAM7DYR5wDFwc++bLc9aOIVUdYyqXljeFubnpeIW5R2EWxi3APhcRFoHlPknrtv1KuBw3C/JTBFpHlBmMm7x+HOBo4G9gKlB65a+BvQHRgIn+J6/EmY9TTVZ6o7Ga9u2bXzyySdceeWVJQI1v1atWlX6nE8++SQbNmzgtttui0QVjamsqcDtItI3cKeI9AFuxSWONxE0a9asokQp942/z/Ks1eaHqeqIwNci8jdgJ3AULuGg4HK53a+q7/rKXIAL2P4KPCMiLYCLgQtVdWbAeVYBxwGfikhvXIA2WFW/8pX5O/CFiOyvqr/W+MUaUwPGj4fFi2v3M/v1g8mTwy+/YsUKVJXevXtXWPbWW29lwoQJJfbl5+fTsWPHotc//fQTd911FwsWLMDj8WBMFNwBDAe+E5FvgTVAZ2AgsBKXON5E0MyZM13KfSDfm2951sIpJCLnV1RGVf9bhc9vjmvd2+573Q3oAHwWcN5sEZmHa417BjgMiAsqs1pElvnKfAqkAHuArwI+az5u0flBQKlgTUQuw3W9Eh8fX4VLMYGsG7Txqkyr6nXXXcfFF19cYt8LL7zA66+/DkBubi7nnHMODz/8MN26dYtoPY0Jl6puEZHDgetwQVs/3JrZ/8Ilht8Zxeo1SN26dQMvXLIQztqqdN+xI9pViqpwW9amlLE/8H/lqgRrjwKLAX/7Zgff48agchtxf8X4yxTiflGCy3QIKLNZA741VFV9EyM6EIKqPgs8C5CYmGh9eNVkEwxqRmVauKKlV69eiAjLli3jtNNOK7dsmzZt6NmzZ6l9fuvXr+fnn3/mwgsv5MIL3WgLr9eLqhIbG8v06dM5/vjjI38RxgTxjXe+w7eZGtayZUsu2QbPTgVQePBB6NEDLrss2lWLinBTd3QLsQ0A7gKWA0dU9oNFZBIwGDhDVQuDDgd/00uIfaVOGVQmVPlwzmOMqYbWrVszYsQIHn/8cfbs2VPq+I5K/IXcuXNnfvrpJxYvXly0jR07lp49e7J48WIGDRoUwZobY+qKtWvXMi7LfWn7+2l2PTU5ijWKrnDXBl0VYvcqYJFvnNl1uDFlYRGRfwPnAMeo6h8Bhzb4HjsAqwP2t6e4tW0D4AHaUnKdtva4uSP+Mu1FRPyta756tqN0q50xJsKefPJJBg0axIABA7jnnnvo27cvqsqcOXOYOHFi2Kk74uLiSuVUa9++PQkJCZZrzdQqETkYN156f6BJ0GFVVcu1FkFr1qyhWdBomsw929krOtWJunBb1srzBXBiuIVF5FFcYDdMVX8JOrwSF2gNDyjfBDfj0z/+7DsgP6hMF6B3QJl0IAk3ds0vBUik5Dg2E2FiSxcY3HiTRYsWMXz4cG688Ub69u3LsGHD+Oijj3jmGUv2buoXETkCt5D7SGAELq1Ud1yGg57Ymi0Rt2bNGloWlLytLbOjVJk6QKqbYkFEbsHlQescRtkngL8BpwI/Bxzao6p7fGVuxE2FHgP8hptlMwTYX1V3+8o8hcufdgEuv80k3C/PYf4uVRGZAXQBLsX9Ij0LZKjqyRXVMzExUTMzMysqZkIY98w4ntrwFJ+d+BnDBwyv+A3GmEZHRLJUNTHa9QiXiMzC9cr8DddYMEBVF4nIMFxKqL+p6uzKnNO+Z8p34IEHMnPFr3TO9wJu/JIceCAsXVr+G+uZcH8Xwp0NGmpAZTxwMK5V7fEw6zXO9zgraP9dwATf8wdxiW6foDgp7vH+QM3nWlyOtjcpTop7ftDYt9G4pLj+WaMfAVeGWU9TRf7g3yYYGGMakL64xgH/f2weAFWdLSL3AhOpwthtE1p6ejrLli3jd4+bWaj4mi5POim6FYuicGeDTgixLxc3bu1fuB/UCqlqhU3FvjFmE8r4TH+ZHFzS3KvKKbMNOC+cepnI8XeDWuoOY0wDEgdkqqpXRLYBHQOO/YpruDARMm3aNI4EBvmaX/wzA6Vly+hVKsrCnWAQibFtphHwt6hZy5oxpgH5neL0UT8CF4nIVN/rCymeHGcioGXLlqRSPBDQ/23yu+ygR3SqFHVhBWEiMkRESq8b444liciQyFbL1FdFLWs23tYY03BMxU0mALgPN9FgFy6h+19x46ZNhOzcuZN5IkXRmuDWdF+9cnEUaxVd4XaDzsHNpvwmxLH9fcdtHRhja4IaYxocVb0z4PnnInIkbn3qZsAnqvpZmW82lfb111+T2bcvuwr20Grp7wDkxUKbkWdEuWbRE26wVl4zSQJFK3gZ42MNa8aYBkBE4oBRwI+quhJAVb8Hvo9qxRqo+fPnM2/ePK4ZOJCWXy0p2r/8zGH0PbVxrl4A5QRrIpKMyyPjNyBEV2hT4CIgvAyXxhhjTD2iqvki8hZwAi4XqKkh6enpHHfcceTm5hI3fz54vUXHYmfN4acPnqVPIw3YymtZuwC4Eze2T4H/ULK9xD+btgC4oqYqaIwxxkTZH7hVckwNSktLIzc31z0HVEDVBRoHbFYKz7kS5vSBlJTyTtMglResTcHdLwFm4wKyn4PK5AK/+dJkGGMrrxpjGqIHgVtFZLaqbq6wtKmSoUOHFj1flJDA7iYeMvP30HGPbzZkfgGkpVmwFsi3HugqABE5BvjOv8qAMRWxPGvGmAZkGNAaWCkiC4D1lPzTVFX1gqjUrAEpLCxEVTn99NO54YYbaJJ6NBktoU02xBf6bnhqanQrGSVhpe5Q1bkWqJmwWMua8Vm7di2XXXYZXbp0IT4+ns6dO3PppZeyZs2aojITJ07k8MMPZ6+99qJdu3acfPLJLFmypJyzGhMVg3HLTG0GevheHx20mWp68sknadmyJa+88gop2dkk5BVy8Cb3tfJzW1+2gXffhfT0aFe11oU7GxQRGQGMxaXqaBJ0WFW1seaqM6FYw1qjtnLlSgYNGkS3bt14+eWX6dWrF7///ju33norhx9+OOnp6SQnJ5OWlsa4ceM4/PDDUVXuuOMOjjvuOH7++Wdat24d7cswBgBV7Rbpc+bnt+HhhyE+Pvzt55/h++9h2DA49lhISICG0okxbdo03nrrLc4880yaNWsGU13OYQ9unsHCTnDgFuCRR+Df/4ZzzoFTT4UBAyA5ueHciDKEtZC7iIwCPgY+B4YDn+DyyxyF6yr9QlUvrMF61ipbYLfqLn/icp7e8jQzT5nJcf2Oi3Z1TJSMGjWKH374geXLl7v/eH2ysrLo1asX/fr1Y9q0aaXet2fPHlq0aMEHH3zAySefXJtVNrWoPizkLiLzgMtU9ZeAfcOAr1W12l8QIgMUFlbrHHFxsNde0KKFeyzveXnHEhJCnz893Q0RS02t2WFi6enpDBkyhIKCApo0acLs2bNJmT0bbruNAoE8D7zcDy77DjxFq7N6oNCXNaxtWxe0DRgAhx/uHjt1qt2LqKKILuQO3I5bWP1aXFPwbaq6SET2Az4FZlS5pqZBsjFrjde2bdv45JNPuPfee0sEagDNmjVj3Lhx3H777Wzfvp1WrVqVOL579268Xm+p/cZEwWBgL/8LEfEAM4HDgUXVPXnTpr+wcSPk5YW3vfaa27xeiImB446D/v1h507Ytat4W726+PnOnZCfX3FdEhJKB3EFBS7O8XohNhbuvRdGjYJu3SAxwmH2lClTKCgoACA/P5+0tDRS8vIoBO4aCp/7+u0u+jEGT6G4ZsYZMyApCb791m0LF8LEicUBXKdO0L07LFjgLiIuDt54A0aMgKZNI3sBtSDcYO0A4A7cig/qf5+q/iYiE3DB3Fs1UUFTv9gKBjVr/CfjWbxhca1+Zr8O/Zh8wuSwyy9fvhxVpXfv3iGPH3jggagqy5cvZ+DAgSWOXXPNNfTr14+UOvgXsDFEcICHiNK8efjlW7Z0w7Xy8lysMmFCxQ1FqpCbWxy4BQZxFT1fsaI47snPhxtvdBtAu3YuaOvWzcVD/ufdusE++7i4KNwGrTVr1vD2228jIsTExBAfH09qaircdBNZezXh8x45LOjqyt5562Dujzuh5EkPOwzGjnXPs7Jg8eLiAO7TT13UCe5GnHaae96smWuNa9fOPZa3tWsHrVu7i4qicIM1L1Cgqioim4F9KF56ah002rVVTRBbwN34ldW66g/og49fd911fPnll3z55Zd4PLZ6nTGBUlJg1qzK9eiJQJMmbmtfySxx6eluXFxenotT/vMf15C1cqXb/vjDxUPvvlscD4Fr9WvXDjZvdsFifDx89hkMCbGCeE5ODqeddhoFBQW8+uqrZGRkkJqaSgrAF1+QpMqsl+HYC2BBV5gc8zWnnHc/KV3LuPhmzWDQILcFXkRurmsevPZa13S4ZUvJbfly97hrV9k3pGXL8gO64H0tW7qbEXhDq9EdG26w9iuQ7Hu+EBgvIvNxCXGvBzIq/cmmYfLHatYLWiMq08IVLb169UJEWLp0Kaeeemqp48uWLUNE6NGj+G+8a6+9ljfeeIM5c+bQvXv3Uu8xJkpC/fUZtb9IU1Jqb9hVuMFhQQGsXVscxK1cCR99BBs3uuO5uTB8OJx4Iowc6bYuXeDZZ3/krrvmsW5dLB988AqnnHJK8Umvvx5UESDBK6RmKAu6Ql5hHhPSJjAhdULZAVtVLsIvL690IBdqW7vWteBt3uwuMJSYGGjTxgVu8fGwZIlrqvR44IQTYO+9SwZzFQg3WHsV8Pdp3ImbaOCff18I/DXsTzQNm++/MRuz1ni1bt2aESNG8OSTT3LttdeWmmDwxBNPMHLkyKLZntdccw1vvPEGaWlpHHDAAdGqtqmDRGQc8A+gI7AUGK+qX5RRtgnwNNAf9301X1VTQ5QbCkwCDsL1DD2oqk+XUYW7RGSL/62+x3tEJDgRfIPMsxZOcBgbC/vu6zZ/CrSRI4tb5Twe93rhQnj/fXe8U6ds1q3rDRwIXMTGjb+XPOnu3e7R44G4WNJ7KJCHonz2x2fMXTWXj879iON7HB+Zi/CLj3dj3fyTEyqi6rpeywvsNm+GRYuK+5QLC+GLL6B5c/f+MIUVrKnqEwHPvxORPrh10poBn6tq8MoGppHyd4NasNa4Pf744wwaNIjjjjuOe++9t0TqDlXl8ccfB+CKK67glVde4YMPPqBVq1Zs2LABgKSkJJKSgpciNo2JiJwNPAqMA770Pc4QkQNVNdR61B4gB3gct/B6yxDn7AZMB14EzsNNInhSRDar6rtBxf+kuJHCbxUuyAtm4z8ChGrQUoVly+C551bz+OOFwL64+DeGN97YxWX+JT/T0+HDD6FnT7joIjypqUzsArfNuY3ZK2cDkFuYy4j/jeDAdgcyqMsgUrqmMKjrIPZrsx8xEn5rVbWJuNkWiYkuWi1LYJ9yfDx88klxABnmd2WFqTtEJB64HJilqo0iW6Wl7qi6Sx+9lOd3PM/sM2ZzzMHHRLs6JopWr17N3XffzfTp09m0aRPt2rVj1KhRTJgwgS5dugBlB/V33nknEyZMqMXamtoUTroCEfka+FFVLw3Ytxx4R1VvruC9jwMHB7esicgDwOmq2itg3/PAQapaq7NaGtv3zIYNG7j11lt56aWXSEoazu7dHwAJgNCsWSFjx8bSZ0865085hpi8XDdQbu7coqAmfXU6x/73WPIK84iNieWCQy5gze41pK9OZ3vOdgBaN23NkV2OLArgBnYeyE8bfyItI43U5NTwuk5rShlj1iKWukNV80TkfmBEdeppGgnrBjU+Xbt25bnnniu3jM0eNqH4GgkOAx4OOvQZMKgap07xnSPQp8AFIhKnqmEkujCVkZuby2OPPcY999xDTk4O119/PbfddhtvvrmCd9/dypFHdmLatP2YNAluYg5uyXFcuo2AdUBTuqYw6/xZpQIvr3r5betvfLX6K9JXp/PVmq+Yvnw6ABIweNoT4+Gqw6/isE6H0bppa1o3bU2bZm1o3bQ1LRJa4Imp4UlN1RxwGO6YtWVAd2BelT/JNCoWrBljyhErIoEZYZ9V1WcDXrfFdWtuDHrfRqA62bY74MZcB58z1veZ66txbhPgq6++4umnn2b27NmsXbuWk046iUceeYT99tsPgMsu61PU9ZmQ4FZmaOndSQzgJQaNjccTtA5oSteUUq1jMRLDAW0P4IC2B3DRoRcBsD17O1+v/ZoH5z/InIw5ABR4C/j31/8OWVdBaNW0VVEQ17ppa9o0bRPyea0HeT7hBmt3AI+KyHeq+lNNVsjUb0UtJRarGWPKVqCqA8IoF9z0KiH2VVaoc4babyopNzeXr7/+mldeeYUXX3wRr9eLiDBp0iSuvfbaMt93zDHwf573GO+dzM8cwCv8jXQ9homkUJW2qFZNW3FCzxNokdCiqOs03hPPu2e9S4/WPdiWvY2tWVvZlr2taNuaXfx6a9ZWftv6G9uyt7EjZ0eZnxMc5JUK6gJer9u9jiWblnBElyM4rONhxEhMpcbXhRus3QgkAd+LSAbur4/AH2xV1aFhf6ppsKxbyxgTAVtwmQY6BO1vT+nWtsrYUMY5C4Ct1Thvo5STk8OCBQuYO3cuaWlpLFiwgJycnBJlYmJiSu0LlrL8vxxZMAZQupFBGsfwTUFKYC9olZTVdVoZhd5CduTsKBXMhQr0tmRt4betv7E1e2vZQd7XVbuWcIO1QsBmfJqwWTeoMaaqfGOlv8OtRf12wKHhQPCszcpIB04N2jccWGjj1SqWnZ3NggULSEtLY+7cuSxYsIDc3FxEhH79+nH55ZczdOhQmjRpwmmnnUZeXl7xigShqMIzz8CVVyK+P/RjySeVNBZ4U8jIcOPyqxuwVWdigSfGQ5tmbWjTrE2l3ucP8rZlb2PSgkk8u/BZvHiJIYYzDzqTk/c7mUJvIWMmjAnrfOGm7kitVC1N42UTDIwxkTEJeEVEvgHmA2OBTrhcaojIRGCgqh7rf4OIHAjE48afJYlIPwBVXewr8jRwpYhMBp4BjgLGAOfW+NXUM4WFhcyYMYP3338fVWXFihV8/fXX5OXlERMTw6GHHsoVV1xBamoqgwcPLrWe76xZs0hLS3MrEoSKtj7/HK67Dn76CQYOhB9/hPx8YmLjyemXCl/Dc8/BK6+4NCD1bQW6wCDv/L7n8/Lil4u6Y8cfMb4ogBzDmLDOF27LmjFhKcqzZoPWjDHVoKpvikgb4DZcUtwlwChVXeUr0pHSSx1OxyXw8vve9yi+c64UkVHAv3EpqdYBV4fIsVaCiLQFmgXmdxORvwMHA5+q6tQqXGJUeL1etmzZwrp161i3bh3r168veh64bdiwAa/XW/S+3r17c/XVVxcFZy1atCj3c1JSUkoHaQUFLh3Hf/7jcqmBS9ExaZLL5p+Whic1lXZzUuBr1/CWl0e1u0OjLRLdsWEHayLSGbe01BCgDXCyqi4RkfFAuqpWsSfWNCRlrftojDGVpapPAk+WcWxMiH3JYZxzLm6Vg8p4EbdqzzgAEbkduAvYDowTkb+q6puVPGelpaenl9lapaps3bq1zODLv3/9+vUUBC7m6dO2bVs6duxIp06d6NOnDytXrmTu3LmoKh6Ph7/97W/cfHO56e2CK+uirMGDXcT19tvw3nsuo3/gouheL8ybBzffXBSRHYNbGaGgwOWMbVO5Hsg6qbrdsWEFayJyEPAFbuxaOnAorqkZ3F8xA7ElpwwUTzupxSTSxhhTwwYALwe8Hgvcp6q3ichjwHVAtYM1r9dLTk4OWVlZZGVlkZmZWfR84cKF3HTTTeTn5xMbG8spp5yCqpYIxvLy8kqds3Xr1kVB2AEHHECnTp3o1KlT0b5OnTrRoUMHEhISSrwvPT2dY489tuJxZ8Fycty6UmPGQH5+8ZJKzZrBySfD//0ftG7tFgv1Z/QPTtGRAldf7RrcCgth/Hjo06d+t65VV7gta4/gcq2NwC3nEfgT8RXwQITrZeopf8tarS75YYwxNas1vlmoInIwbkapP3j7ADi/sifMycmhd+/eJQKz7OzssN6bn5/P1KlT6datG506dWLIkCEhg7COHTvSpEmTylYNcN2YX0+ezNZ336XNGWfQxx8pqcK2bfDHH/D778Wb//XatSXXvBSBc891A9AC1gmuaIH1li2LP64hdIVWV7jB2mDgXFXdIyLBGeA2UnoqtGmkvOrGOMTEWLBmjGkwtgJdfM+HAetUdbnvdRxV6EsQEfr27UuzZs2KtsTExBKvA/f9/vvvXH311RQUFBAfH8+sWbNCD9yvCq8X9uxxC6j7twUL6HPDDa51bPZseOMN2LnTBWQ7d5Z8f8eO0L07DBsGPXq45rAHHnD9mPHxcOWVJQM1qDCj/3HHwR13uOceT6nGt0Yn3GDNW86xtkB4fw6YBk+9NmbNOBs3buS+++5j6tSprFmzhrZt29K3b1+uuuoqRo0aBcDXX3/Nv/71L7788ksyMzPp1q0b5557LjfeeGOVWwSMqQGfAxN8Ew2ux7Wm+R2AW+C9UhISEnjzzfB7To899lhSoLilq39/2LKlZIC1ezfs2lV6X0X7K1qjtKDAzdo8/HAXYPXo4bbu3d0WHIgBnHBCuS1n4RBxLWuWvjP8YO0b4ELg4xDHzsJNqzamuBvUWtYatYyMDI466iiaN2/OxIkTOeSQQ/B6vcyaNYuxY8fy559/8tFHH3HmmWcyevRoPv/8c9q0acNXX33FDTfcwKxZs/j888+Jj4+v+MOMqXn/BP4HTAS+xU0u8BsNfFnZE3bOzYUzz3QtV/6toKDk68B9e/bQZ6MvH/Bnn8Hf/x7eBzVpAs2bw157ucfmzWHvvaFnz5L7/Jt/36pVLrWGv3Xso48qF3RVcy3M//63OEjLz3evrRu0YvcAn4vIZ8BruGHkx4nINcBpuBmixlg3qAFg3LhxqCoLFy4kKSmpaH/v3r0ZPXo0WVlZXHzxxYwaNYqXXnqp6Pi+++7L/vvvz4ABA3j00Uf5xz/+EY3qG1OCqm7EJc8N5TjcWO5KaeL1wrJlbmakf4uNdQtlJiWV3BcXB7/8Aps2uQhGxPUTnnRS+QGX/zxVdeih1W4di5Tnn3fV2bq1ctXxT0qN9iVUtx4S7vJAInIiMJmSeW0ygCtUdUblP7ruSkxM1MyKmoVNSOfecy5veN/gu0u+o3/nys6ONw3Btm3baNu2Lffeey+33HJLyDLvv/8+p59+OvPnz2fQoEGljg8fPpwtW7bw/fffh3i3qe9EJEtVE6Ndj+rwJeDtjUtdta6y76/090x6Ohx7bPEMyvqYKbYS0tMhxH8NgItfTz4ZOlQwWn7DBvj4YzeEzuMJ7z01obx6PPVUeL8LYedZU9VpwDQR6YlbS22rqv5axbqbBsq6QWvY+PGweHHtfma/fjB5ctjFV6xYgarSu3fvMsv89ttvAGWWOfDAA3nuuecqU0tjaoyIPA7EqupY3+vTcak6PMAuERmuqt/WaCVSUiqcQdmQlHd5BQWuJzjUULlAWVmubGXeUxMiUY9Kr2CgqiuAFZV9n2kcLHWHCbe1vqJz2CQVU4eMpOQ4tbuAqcAduNRWdwIn1XgtqjkOrL7p0MG1SgWKiXE9xTNnVnwrghsjw3lPTSivHuH+N1eZFQx64Zb9SAE6A2txOdbu9QVwxqCFXrrvtNmgNaYSLVzR0qtXL0SEZcuWcdppp4Uss99++wHw888/c9RRR5U6vmzZMnr16lWj9TSmEjrghv0gIl2Ag4CLVfUnX1LcF6JYtwZr/XqXFWTDBjf87pFHKjdmra40RkaiHmGNWRORVNyaa9nANFxutb2BE4FmwAm+JTwaBBuzVnXPntify6Z/z+9TX6HHiedFuzomSkaOHMkPP/zAb7/9VmKCAcCOHTuIi4tj3333ZfDgwXzwwQclji9atIgBAwbwwAMP2ASDBqq+jVkTkS3A+ao6XURGA08ArVXV6/9+VNVKdWzZ94yB8H8Xwu2regS3IO6+qnq+qv5DVc8HkoHFvuPG0Gv1FgDif8+IbkVMVD355JOoKgMGDODtt9/m119/5ZdffuGpp56ib9++JCYm8txzzzFt2jQuuugivv/+e/7880/eeOMN/vKXvzB48GCuueaaaF+GMX6LgCt8qxdcAcxUVX/+0W7A+qjVzDQK4XaDHgicrap7Aneq6m4ReQB4PeI1M/VSJMYrmfqvW7duLFq0iPvuu48bb7yRtWvX0qZNGw455BCeeeYZAE477TTmzZvHv/71L4YNG0ZWVhbJyclccskl3HTTTZZjzdQltwKfAD8AO3Brg/qdistFWilZWVkqIg0xoXwsUHql+IYj0tfXNNwPDccaihduDxaPG79mTHGmaQvaGr2OHTvyn//8h//85z9llklJSWHq1Km1WCtjKk9VvxWRfXCrFSxX1V0Bh58Flod+Z7nnbJCzsERkoaoOiHY9akq0ri/cYO0B4C4RSVfVosBMRDrjZsHcVxOVM/VPUcuaTTAwxjQgqpoJfBdi/7QoVMc0MuEGa0OB5sDvIrKA4gkGR/qep/oGWQKoql4Q4XqaeqIoWLOWNWNMAyMihwD7A6UWrlXV/9Z+jUxjEW6wNhgoxA2i3Ne3QfGgyqMDytq3dCOmWLBmjGlYRKQlLhNCCu47zt91EPgfnQVrzrPRrkANi8r1hRWsqWq3mq6IaRisG9QY0wDdB7TBNUx8gVsTeydwES6AOyd6VatbVLVBB2vRur4GOcDRRI96rUXNGNPgjMAFbAt8r9eoapovhdXngOWZMTWqUstNiUhXoCuh++tnR6pSpv6yblBjTAPUEfhDVQtFJAc3htvvPeCN6FTLNBZhBWsi0h14FRjo3+V79PfdK25BW9PIWcOaMaYB2gC09D1fhev6TPO97hmF+phGJtxu0OeBfYDxwAnAMb5tWMCjMajXl9TbWtaMMQ3Hl7gADeAV4E4ReUZEngAeAj6NWs0iQERai8j7IpIpIqtE5K/llL1WRDaIyE4ReVFEEgKOXSkiC0UkV0SmBL0vXkTeEZEMEdGADBL+4xNEJF9E9gRs3evR9R0pIjNFZJuIbBaRt0WkY8BxEZEHRGSrb3tQKrGIdrjdoIcDY1T13XBPbBonb9EKLMYY02DcBXTyPX8IN9ngbNza2B8BV0WpXpHyBJCHS8nVD5gmIj+o6tLAQiIyArgJ10CzDngfd29u8hVZB9yLG+MXKjP/l8Bk4O0y6vGmqtbEotK1cX2tcDNFP8WtcPA48BKugQvgMtxqF4fgeiNnAn8AT4dzAeG2rK3BXagx5fJ6LVgzxjQsqvq7qn7he56vqterahdVba2qf1XVrdGuY1WJSCJwBnC7qu5R1S9xAejfQhS/AHhBVZeq6nbgHmCM/6CqvqeqHwCl7oeq5qnqZN/5CyN/JaHV4vXNUNW3VXWXqmbhgrWjgs79iKqu8S0u8EjguSsSbrB2H3Cj76KNKZMFa8YYU6/sBxSq6m8B+34ADgpR9iDfscBye4tImwjV5WRfN+JSEbk8QueM1vUNAQJb7kKdO1QdQgo3z9orInIAkOFbwWB76SK2aoGBQgvWjDENgIhcVJnyqvpiTdWlhiXhcsYF2knJGa9llfU/b06I1qZKegvXjbgROAJ4V0R2qOrr1TxvrV+fiPQF7gBOqeDcSSIiqhUP8g53NugY4GZc02V/SneJ2mhyA4DXa6k7jDENwvMUf7dVNBBcgfoarO0B9gratxewO4yy/uehylaKqv4c8PIrEXkUOBOobrBWq9cnIj2BGcA1/q7zcs69J5xADcKfYHAXbqDdxaq6I8z3mEao0DfBQG2igTGm/tsDvIObAboyynWpKb8BsSLSS1WX+/YdQskuPL+lvmNvBZTbWENj9gKX9aqOWrs+EdkXlyT5HlV9pYxzf1NBHUIKd8xaG+BJC9RMRbyFFqQZYxqEbsDDuCWmPset/XkssE1VVwVv0axodahqJi6x790ikigiR+G674KDDXD34GIROVBEWgG3AVP8B0UkVkSa4PKuekSkiYjEBhxP8B0HiPcdF9+xU0SklS/FxUDgauDD+nJ9ItIZmA08oaqhZnj+F7hORDqLSCfg+sBzVyTcYO1LoHe4Jy2PiAwRkY9EZK0v18qYoONTfPsDtwVBZRJE5D8issWXN+UjEekSVKaViLziy5Wy0/e8ZSSuwZStKHWH9YIaY+oxXxB2j6ruhxssvgyXtmODiLwuIiNFpKEs2TgOl4piE67b8XJVXSoi+/jyne0DoKqfAA8Cc3DJgVcBdwac5zYgG5fq4jzf89sCjv/q29cZl+IiG9jXd+wcYAWuy/G/wAOq+nI9ur5LgO64HHxFueIC3vsM8DHwE7AEmObbFxYJp7tURPbHNQs+CHxC6QkGaJj9XiIyChgMLML9g4xT1SkBx6fg/iEDp9Xmqeq2gDJP4SLjC3CD/ibhsksfpqqFvjIzcIl8L8WFDs/jlgs5uaI6JiYmamZmZjiXY4K80DqGi7crqx+8ja7/uCfa1THG1EEikqWq9S67gIjEAyfjvntGAm+rapkJVo2JlHDHrC3zPf63jOMa7rlUdTowHYoCs1ByVXVDqAMi0gK4GLhQVWf69v0NFwEfB3wqIr1xiegGq+pXvjJ/B74Qkf1V9ddw6moqJycnp2iCgUZipIExxtQtbYBkXGuQB9gS1dqYRiPcYO1uardja7CIbAJ2AHOBW1V1k+/YYUAc8Jm/sKquFpFlwCBc02oKbmDoVwHnnA9k+spYsFYDdu4MmJVss0GNMQ2AiDQFTsf19hyHSxL/KnCW/eFvaku4rWETargegT7BDQZcifsL5l5gtogcpqq5QAdcCpHgv2g2+o7he9wcOCVWVdUXAHYgBBG5DLccBPHx8RG7mMZkx44dBNzwaFbFGGOqRUSOwwVop+EaK94DhqvqnKhWzDRK4basFRGRJFxT8DpVzY90hVT1jYCXP4nId7guzhNxvyxlVo2SrX+hooXgMoGf+ywuIR+JiYkWaVTBzp078Q9cFK/dQmNMvfYZsAuXuuM9IAu3HvewUIVVdXYt1s00MmEHayJyEq479BDfrsOBRSLyPDBbVV+rgfqhqutEZA3Qy7drA26sQFtgc0DR9sC8gDLtAzMD+6YHt8O1wJUrK2s/Wrd2z0WKt5gY8HhKPgbvi411W1xc8fPg16GOxcdDkybhbQkJxc+bNoWWLaFFC3euaNqxYwcF/rFqBbW29JsxxtSUvXDrNwau0BM4ItefC0xx30vG1IhwVzA4FXgXmAXciJsV6rcS94NcI8GaiLTFzQ5d79v1HZAPDPd/pi9tR2+Kx6il45Z2SAnYlwIkUnIcW0ixsTs57zzXkxe4eb1uKyws/bywsHjLz4eCArfl50N2NuzeXfw68Jj/MS8PcnLcVtUexKQkF7iFu3XuDN26ueAvEgKDNSksiMxJjTEmOo6JdgWM8Qu3LeZO4CVVvcSXAC4wWFuCy2ESFl83ak/fyxhgHxHpB2zzbRNwgeF63Ji1ibjcKO8DqOpOEXkBeMg3Bs2fuuNHXOJCVHWZiHwCPCMil+L+8nkGmBrOgND4+I089li4VxRZqi6A8wduFW1ZWbBzJ+zYUXLbvh3WrIElS9zrnTtDB4Ei0KUL9Ojhtp49i5/36OFa7MK1Y8cO/O1pYslxjTHVJCLjgH8AHXHZ3scHLeETXL4P8DgwEPd98gwum7y/hyUVl0MrWG9V/SVwh6rOjcAlGBMR4QZrvYF/+p4Hf+Vvx41hC9cASv6y3OXbXgYuB/oA5+Pypq33lT1LVQPX5roWKADexCW6mwWc78+x5jMaeIziWaMfAVdWop5RIeK6RePioHmoZWaryOt1rXuBwdzq1fD778XbRx/B5s0l39epE5x0Epx6KgwbVn4r3KpVq2jhTxFZaN2gxpiqE5GzgUdxjQFf+h5niMiBqvpniPJ7ATNxw2EOB/bHZYjPBB4JKn4QLpjzC/qfz5i6JdxgbRdujFgoyVTiB11V0yh/va8RYZwjB7jKt5VVZhsuw7DBjalr0cJt++5bdrndu0sGcN9+C6+9Bs8+64LHkSPhtNPcY3Cr2/LlyxmcmAhbMxEL1owx1XMdMEVVn/O9vkpETsD9UX9ziPKjgWbABaqaDSzx5dy8TkQmBS2YvUlVLUeaqTfCXSpjJnBz0HJNKiIJuNaqGZGumImO5s2hXz844wz45z/h7bdda9u0aXDOOZCWBueeCx06wNNPl+xaXbFiBU39zYEFNmbNGFOmWBFZGLBdFnjQt1LAYQTk0/T5DJcrM5QU4AtfoOb3KdAJ16gQaKGIrBeRWSJiY9NMnVdmsCYif4iIf+bnrbj8ZL/ilm1S3NpYi4EuuHFmpoFq0gRGjXKta+vWwZdfwtChcPnlMHq0a41TVZYvX05CS9fcFpOVXcFZjTGNWIGqDgjYng063hY3uzJ49n5gPs1gHcoo7z8GbmjN5cAZuES3vwKzRGRIFa7BmFpTXstaMpAAoKoZQH9gKm4WZiFuYdsFwBGquq5Ga2nqDI8HjjoKpk+Hf/0L3nwTBgyAr7/ewp49e5CO7QGIX7U2yjU1xjQAwWOky8yVWU75ov2q+quqPq2q36lquqqOwyVivyEitTWmhoSdmUtV1+DW5DSGmBi45RYYNMiNYbvqKvej1LJlSwCarCw1/tcYY8K1BdcoENyK1p6yc2VuKKM85bwH4GvgnMpW0JjaVNGYNUtDb8qVmgp33gkLF7YCRtCujZuHkrBydVTrZYypv1Q1D5dTc3jQoeGUnSszHThaRJoElV8HZJTzcf0ozuNpTJ1UUcvaXSISzowZVdULKi5mGqLLL1f++c8/iY9/kmZNJwHgycmF3NzIZdw1xjQ2k4BXROQbYD4wFjdZ4GkAEZkIDFTVY33lX8PlBJ0iIvcC++HGVt8VkGdtPC5wWwrE4zIGnIobw2ZMnVVRsNYPyA3jPNYC14gtXfo9+fmPkJ//KhvXBWRlWbECDjooehUzxtRbqvqmiLQBbsMlxV0CjFLVVb4iHYEeAeV3ishw4AlgIS4H6CO4oM8vHngYtypONi5oO1FVp9fw5RhTLRUFa6eq6je1UhNTb7366qvExn5OQQGsD5xq8t57FqwZY6pMVZ8Enizj2JgQ+37CTX4r63wPUnIFHmPqhXDzrBkTUmZmJq+//jqjRh3JwQfDhrWuZW3PgT3hhRfc0gnGGGOMqTIL1kyV5efnc9ZZZ7Fx40auueYajj4aNvha1tadexKsWgUzZ0a3ksYYY0w9Z8GaqRJV5bLLLmP69Ok89dRTDBs2jAEDoDDXtaytP+5IaNsWnnuugjMZY4wxpjxlBmuqGmPj1UxZbr/9dqZMmcKdd97JZZe5lWK6dAGP1/1IZccBF1wAH34Izz8Pc+e6leNtzVBjjDGmUsJOimuM1+tl2rRpTJ48mdmzZ3PppZdy5513Fh3v0gU8hS5Yy/LmufWoXn4ZLr20+CTx8ZCcDN27Q9eubkkEVTe2TbX+Pg91rL4osb51HVZf6mmMMRFmwZqp0O7du3nppZd47LHH+P333+nSpQsPPfQQ48ePR6Q4VUdgy1qONw969HCLif75J/zxB6xc6R792/ffuy9gEbfFxIR+Xt6x6j73eKp/nvKO1Rf1pa71pZ6mfFOmRLsGxtQrFqyZkDZv3syPP/7ItGnTeOGFF9i1axeDBg3ivvvu47TTTiMuLq7Ue/baCxJi3P5t+Tvdzrg4F7T16FGqvDGmkbJgzZhKsWCtkcvLy2PZsmX8+OOPJbYNGzYAEBsby1lnncU111zDwIEDKzxf81i30ssfOzNqstrGGGNMo2HBWgOkquzevZstW7awdetWtmzZUur5li1b+O2331i2bBkFBQUAJCQkcNBBB3HCCSfQt29f+vbtS79+/WjTpk3Yn90s3o3V+m3H7zVybcYYY0xjY8FalBQWFlJQUEB+fj4FBQXk5OSQnZ1dqS0rK4vMzEy2bt1aKijLz88P+bkej4fWrVvTpk0bevbsycknn1wUmPXq1YvY2Or9SDSNL8QL/L7jj2qdxxhjjDGOBWsh5OTkMHz4cHxr/xY9er3ekFthYWHR8/z8/KIArLxHrebMtpiYGJo2bUqzZs1o27Ytbdu2pWfPnhx55JG0bduWNm3aFO0PfN2iRQtiYmouvV7T+EIKJYaV21fiVS8xYqn8jDHGmOqwYK0MWVlZAEWzHUUEEcHj8RAbG0tMTEyJzePxICLExcUVbbGxsSUey3oeGxtLQkICTZs2DXuLi4srMROzrmgWl0++eMgtzGXtrrV0bdE12lUyxhhj6jUL1kJo0qQJ8+fPj3Y16qUEySNfEoB85q2ax+i+o6NdJWOMMaZesz4qE1HxMfkUaBP2a7Mfj33zWLSrY4wxxtR7FqyZiIqXfPI1jqsHXs03a7/h6zVfR7tKxhhjTL1mwZqJqHjyySOOC/pdwF4Je/Ho149Gu0rGGGNMvWbBmomoOPLJ0ziS4pO4+NCLefvnt3nx+xdZuX1ltWfAGmOMMY2R2BdoaYmJiZqZmRntatRLSw88E5Yt40DvUlbv+pOjXjyKNbvWANBlry4M2XcIQ/YZwpB9h3BA2wPq5IxWY0zNEpEsVU2Mdj2MqS8sWAvBgrWqW9FrJNtWbKV//jfExoJXvSzdtJQv/vyCeavmMW/VPNbvWQ9Ai4QWtGraiqaxTWkW14xmcc1oGhfwPLaM52WVCdrfNK6p5Xkzpg6yYM2YyrFgLQQL1qpudbchrMjwkJI9hyZNSh9XVX7f/jtfrPqChesWsid/D1n5WWTnZ5OVn1W0ZRcUv87Ozya3MLdK9WkS26RawWC4gWKcp/TC9saY0CxYM6ZyLM+aiai4giyy2BvfcqOliAg9W/ekZ+ueXHjohWGft9BbSHZBdomgLjigq9T+gmx25OwIGSgqlf8DxiOeEgFd09jQrXpldfsKpfeHKhuqXGXK1sTnN8RrqqnPN8aYqrBgzURUfF4mmSSWGaxVlSfGQ1J8EknxSZE9cRBVJa8wL+ygr6xj2QXZpSZUlBUEhmrdDlW2rFbwcMvWxOc3xGvyly3171fNzzfGmKqyYM1EVGx+Flk0i3iwVltEhITYBBJiE2hFq2hXx5gGSS62VkdjKsNGX5uIis/dzR6S8HqjXRNjjDGmYbBgzUROXh5Nsrazkb0tWDPGGGMixII1EzmbNgFYsGaMMcZEkAVrJnI2bHAPdKCwMMp1McYYYxoIC9ZM5FiwZowxxkScBWsmcta7lQk2sjeWucAYY4yJDAvWTOR8/z15TZqzhi4WrBljjDERYsGaiZz589nS80i8eCxYM8YYYyLEgjUTGbt2wU8/sbnXUQA2G9QYY4yJEAvWTGTMng2qbNlvEIC1rBljjDERYsGaqb68PLj5ZujenS0HDgGsZc0YY4yJFFsb1FTf5Mnwyy8wdSq6OwGwljVjjDEmUqxlzVRdXh488ADccQf85S9w4omIb31mC9aMMcaYyLCWNVN5qjBnDlx5JSxbBqecAs8+C2DBmjHGGBNh1rJmKrZ7twvOJk6EU0+FTp3g2GMhJwemToUPPoD27QEL1owxkSMi40RkpYjkiMh3InJ0BeX7iMhcEckWkbUicoeI/3+lojJDfefKEZE/RGRszV6FMdVnLWuNldcLO3fCli2webN79G+BrzMyYOnS4uhrv/1g+HA46ig4/3xo2rTEaS1YM8ZEgoicDTwKjAO+9D3OEJEDVfXPEOX3AmYC84DDgf2BKUAm8IivTDdgOvAicB4wGHhSRDar6rs1fU3GVJUFayG0ys+Hxx8vfSA4Agn12r/P/7ysfaHKFBa6IMrrLX4e/Bi8r7AQcnNLbjk5pfcFHy8oKPsGJCRAu3bQti3ssw+ceSYccQQMHAitW5d77/zB2hNPwMEHl1vUGGPKcx0wRVWf872+SkROAC4Hbg5RfjTQDLhAVbOBJSLSG7hORCapqgJjgXWqepXvPctE5AjgBsCCNVNnWbAWQof8fLjqqooL1hQRiIkBjyf0Y/C+hITirUkT99iqVej9gVuLFsVBWeDWrFlx1FVJnTq5tz7zTITviTGmIYkVkYUBr59V1Wf9L0QkHjgMeDjofZ8Bg8o4ZwrwhS9Q8/sUuAdIBlb6ynwW9L5PgQtEJE5V8yt7IcbUBgvWQvitaVP4s1QruxMcxIR67d/nf17WvsDX/iAsJqbKgVJdkJIC27dDvv2XZ4wpQ7t2FKjqgHKKtAU8wMag/RuB48p4TwdgTYjy/mMrfY+fhygT6/vM9eXX3JjosGAthEIR18JkqqRFi2jXwBjTQASPfpUQ+yoqH7w/nDLG1Ck2G9QYY0xdswUoxLWEBWpP6dY2vw1llCfgPWWVKQC2VqmmxtQCC9aMMcbUKaqaB3wHDA86NBz4qoy3pQNHi0iToPLrgIyAMsHdqMOBhTZezdRlFqwZY4ypiyYBY0TkEhHpLSKPAp2ApwFEZKKIzAoo/xqQBUwRkYNF5HTgJsA/ExTfe7uIyGTfOS8BxlB6IoMxdYqNWTPGGFPnqOqbItIGuA3oCCwBRqnqKl+RjkCPgPI7RWQ48ASwENiOy682KaDMShEZBfwblwJkHXC15VgzdZ1oLWcvFZEhuJw2h+H+SrpQVacEHBfgTuAyoBXwNXCFqi4NKJOA+0voXKApMAsYp6prAsq0Ah4D/uLb9RFwlaruqKiOiYmJmpmZWfWLNMYYUyYRyVLVxGjXw5j6IhrdoEm4v5CuAbJDHP8ncD1wFS4L9SZgpog0DygzGTgDF6wdDewFTBURT0CZ14D+wEjgBN/zVyJ5IcYYY4wxNa3WW9ZKfLjIHuBKf8uar1VtHfC4qv7Lt68pLmC7QVWfEZEWwGZci9yrvjJdgVXASFX91Je1+mdgsKrO95UZDHwBHKCqv5ZXL2tZM8aYmmMta8ZUTl2bYNANN626KMO0Lxv1PIqzVh8GxAWVWQ0sCyiTAuyh5Kyh+bg14srKfm2MMcYYU+fUtQkG/vw3obJWdw4oU4jLwxNcpkNAmc0BM4BQVRWRTZTOsQOAiFyGGyfnf51VlQtowGJxuYhMMbsnpdk9Kc3uSWlNo10BY+qTuhas+VU2a3WoMqHKl3ke37p0zwKIyMIKlkJpdOyelGb3pDS7J6XZPSktaF1QY0wF6lo36AbfY3lZqzfg1owLXg8quEx73xg4oGg8XDvKzn5tjDHGGFPn1LVgbSUu0CrKWu3LRn00xePPvgPyg8p0AXoHlEnHzTpNCTh3CpBI2dmvjTHGGGPqnFrvBhWRJKCn72UMsI+I9AO2qeqfIjIZuFVEfgF+wyVE3INLxeFPfPgC8JBvDNpWXNLDH4HPfWWWicgnwDMicimu+/MZYGpFM0F9no3IxTYsdk9Ks3tSmt2T0uyelGb3xJhKiEZS3FRgTohDL6vqmICkuH+nZFLcJQHnaAI8BPyVkklxVweUaU3ppLhXhpMU1xhjjDGmrohqnjVjjDHGGFO+ujZmzRhjjDHGBLBgzRhjjDGmDmvwwZqIjBORlSKSIyLficjRFZTvIyJzRSRbRNaKyB2BKUB8ZYb6zpUjIn+IyNiavYrIivQ9EZGOIvKaiPwiIoUiMqXGLyLCauCenC4in4nIZhHZLSJfi8hfyjtnXVMD92SoiHwlIlt9ZX4RkRtq/koipyb+PwkoO1hECkRkSajjdVkN/KykioiG2A6o+asxpg5S1Qa7AWfj0nxcikvt8R/czNJ9yii/Fy51yFvAwbjF4ncD1weU6YZbtuo/vnNe6vuMM6J9vVG8J8m4yRxjcKlRpkT7OuvAPXkUuAkYiJv9fCdu5Y2jo329UbwnhwHnAAf5fo/O8/0ujYv29UbrngSUbQX8AXwKLIn2tUb7vgCpuATmB+Lybvo3T7Sv1zbborFFvQI1enFuJulzQfuWAxPLKH85sAtoGrDvNmAtxZMxHgCWB73veSA92tcbrXsSVH4q9S9Yq9F7ElDmG+CRaF9vHbsn7wGvR/t6o31PfPfhTmBCPQzWauL/WX+w1jba12ebbXVha7DdoCISj/tL/rOgQ59R9mLuKcAX6haP9/sU6IRrPfKXCT7np8AAEYmrTp1rWg3ek3qrlu9Jc2B71Wpae2rrnojIob7zza1OfWtDTd4TERmHazW6N1L1rS218LOyUETWi8gsETkmAlU2pl5qsMEabjkqD6EXhQ+5mLtvf6jy/mPllYml9BJYdU1N3ZP6rFbuiYhcAXQBXqlaNWtVjd4TEVkjIrnAQuBJVX26etWtFTVyT0SkD65FbbSqFkamqrWqpn5W1uNa4M4ATgd+BWaJyJDqVtiY+qiuLuQeSZVdFD5U+eD94ZSpy2rintR3NXZPROQMXBLnc1R1VZVrWPtq6p4cjVsO7kjgARFZqar1IYiFCN4TEUkA3gBuUNWVEapftET0Z0XdSjOBq82ki0gycAMwr+rVNKZ+asjB2hbcgO7yFoUPtqGM8lBykfhQZQpwS1/VZTV1T+qzGr0nvkDtFeB8Vf2oelWtNTV6TwICk59EZG/cOK26HqzVxD3piBtA/5KIvOTbHwOIiBQAo1Q1uHuxrqnN/1O+xk1QMabRabDdoKqah1v0fXjQoeGUvZh7OnC0uOWsAsuvAzICyhwX4pwLVTW/OnWuaTV4T+qtmrwnInIW8D9gjKq+E6k617Ra/jmJARKqVtPaU0P3ZC3QB+gXsD0NrPA9L+u8dUYt/6z0w3WPGtP4RHuGQ01uuCnlecAluCnlj+KmlO/rOz4RmBVQvgXur743cFPKT8fNWgqVumOy75yX+D6jPqXuiOg98ZXr59vm4dZh7QccGO3rjeLPyTm4dAbXUDL1QOtoX28U78lVwElAL992sa/M/dG+3mjdkxCfMYH6Nxu0Jn5WxgOn+n5ODvKdQ4HTo329ttkWjS3qFajxC4RxuL/WcnF/AQ4JODYFyAgq38cXcOTg/oq7k9LT7IcCi3znXAmMjfZ11oF7oiG2jJq+lrp6T4C0Mu5JWrSvNYr3ZDywFPfHzk7f79A4ICba1xqtexLi/BOoZ8FaDf2s/BPXwpgNbAO+wHULR/1abbMtGpst5G6MMcYYU4c12DFrxhhjjDENgQVrxhhjjDF1mAVrxhhjjDF1mAVrxhhjjDF1mAVrxhhjjDF1mAVrxhhjjDF1mAVrxlSTiGgYW4aIJPuej6kDdU4Oql9qJd57W8D71tRcLY0xxkDDXhvUmNqSEvT6feAHXIJTv1xc8s8U4PfaqVZY7gWmAT9X4j0vAZ8DtwOH1ESljDHGFLNgzZhqUtUFga9FJBfYErzfJ9S+aPq9jHqWSVXXAmtFZHMN1ckYY0wA6wY1ppaE6gYVkSkiskZEBojIVyKSLSK/isiJvuPX+bpQd4nIhyLSLuicsSJys4j8IiK5IrJORB4JWiS7svUc4avLThHZ46vPHVW+cGOMMdViLWvGRN9ewH+Bh4F1wK3AuyLyBLAfcAWwNzAZeAI4K+C9/wNOBh4AvsItpH0PkAycUdmKiEh34CPgHeBu3ALdvYDulb4qY4wxEWHBmjHR1xwYq6rzAERkHW7M20nAgapa6Nt/MHCViHhUtVBEjgbOBi5Q1f/6zvW5iGwD/ici/VR1cSXr0h+IBy5X1V2+fbOrc3HGGGOqx7pBjYm+TH+g5vOL7/Fzf6AWsD8W6Oh7fQKu5etdX3dorIjEAp/5jg+pQl0WA/nAGyJypoi0r8I5jDHGRJAFa8ZE347AF6qa53u6Paicf79/PFp7XCvYHlyA5d82+Y63qWxFVHUFMAL3f8MrwAYR+VpEhlb2XMYYYyLDukGNqb+2AjnA0WUcX1eVk6rqHGCOiCQAR+HGrk0TkWRV3VKlmhpjjKkyC9aMqb8+AW4EWqjqrEifXFVzgdkikgR8CHQDLFgzxphaZsGaMfWUqqaJyOvAOyIyCfgG8OJmgo4CblTV3ypzThEZixvrNh1YDbQFbsa10i2JXO2NMcaEy4I1Y+q384CrgItwKT9ygQzgU2BjFc73AzASmIgbE7cN+BIYrarZEaivMcaYShJVjXYdjDG1TESSgZXAxbgcb4Ua5n8GIiKAB3gBOFZVu9RUPY0xxthsUGMauxdwM0grM9vzVt97zq+RGhljjCnBWtaMaYREJB7oG7DrV1XdHeZ7OwKdfS/zVPXHSNfPGGNMMQvWjDHGGGPqMOsGNcYYY4ypwyxYM8YYY4ypwyxYM8YYY4ypwyxYM8YYY4ypwyxYM8YYY4ypw/4fd9y3Y0DKTl4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "info={}\n",
    "info['label1'] = {'label':'Temperature','units':' [K]'}\n",
    "info['label2'] = {'label': 'CH4','units':' Mass Fraction'}\n",
    "info['label3'] = {'label': 'O2','units':'  Mass Fraction'}\n",
    "info['label4'] = {'label': 'CO','units':'  Mass Fraction'}\n",
    "info['loc_x'] = 0.45\n",
    "info['loc_y'] = 0.75\n",
    "info['xlim'] = [0,1]\n",
    "\n",
    "info['inset_x1'] = 0.55 + 0.5 \n",
    "info['inset_y1'] = 0.4\n",
    "info['inset_x2'] = 0.3 \n",
    "info['inset_y2'] = 0.4\n",
    "\n",
    "spNumber = 53\n",
    "info['xlim2'] = [IgnTime[spNumber]*0.95,IgnTime[spNumber]*1.05]\n",
    "info['xlim'] = [0,time_simulation[-1,spNumber]]\n",
    "\n",
    "x  = time_simulation[:,spNumber]\n",
    "y1 = solution[:,spNumber, T_indx]\n",
    "y2 = solution[:,spNumber, CH4_indx]\n",
    "y3 = solution[:,spNumber, O2_indx]\n",
    "y4 = solution[:,spNumber, CO_indx]\n",
    "helper.makefigure(x, y1, y2, y3, y4, info, fs_label=16, fs_tick=14)\n",
    "plt.savefig('GRI3_TimeProfile.pdf', bbox_inches = 'tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finalize Kokkos. This deletes the TChem object and the data stored as Kokkos views"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "gkde = stats.gaussian_kde(dataset=IgnTime)\n",
    "x = np.linspace(start=np.min(IgnTime),\n",
    "                stop=np.max(IgnTime), num=N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVYAAAEGCAYAAAA+Ib10AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAsQklEQVR4nO3dd3xUddr38c+V3kEMJRAgVKkWOgqoCMpauXVVBMGyLiuLuvq468q698K6+jyK3nhbUWwLuIINWazsLtKk995bCARCSwgkkHY9f8wBxxAgCTM5M5Pr/XrllZlTv4eBizO/c36/I6qKMcYY3wlzO4AxxoQaK6zGGONjVliNMcbHrLAaY4yPWWE1xhgfi3A7gL8kJydrWlqa2zGMMSFm2bJlB1W19rmWCdnCmpaWxtKlS92OYYwJMSKy63zLWFOAMcb4mBVWY4zxMSusxhjjY1ZYjTHGx6ywGmOMj1lhNcYYH7PCaowxPmaF1RhjfMwKqzHG+FjI9rwypqr07/+jX7c/dWoPv27f+J6dsRpjjI9ZYTXGGB+zwmqMMT5mhdUYY3zMCqsxxviYFVZjjPExK6zGGONjVliNMcbHrIOAMRcgPz+fY8cyKCrKo6SkkMTExkRGJrgdy7jMCqsxFbBv3z5eeuklli1bxrp16zh48ODP5nfvPoY6dToDsH37F6Snf0NiYlOSk6+gTp0uxMae8xl0JkRYYTXmHA4ePMj27dvp0qULAOHh4YwZM+b0/MjISCIjk4mMTEAkgoiImNPzsrM3kZOzhZycLWRkTAegRo0WNGp0Ew0bXk9kZGLVHoypMlZYjSnDihUrGD16NFOmTKF+/fps27aNsLAwateuzcsvv0zr1q257LLLSElJ4fbb55e5jfbtHyMt7Rayszdz4MASDhxYTk7OFtas+V/27p1Fjx6vV/FRmapihdUYL+vWreMvf/kLU6ZMAUBEaN26NYcPHyY5ORmAJ598slzbioxMoFat9tSq1Z6mTe+guLiAffvmsnPnVzRseP3p5QoL85zl43x8NMYtVliNAU6cOMGzzz7L6NGjKS4uJiYmhmHDhvH444/TqFEjn+wjPDyKBg2uo0GD61DV09PXrx/L/v3zueyy31O3bnef7Mu4y263MgYoKiri448/pqSkhGHDhrF9+3bGjBnjs6JamogAUFxcQHb2JvLzs1i48ClWrPh/FBWd8Ms+TdWxM1ZTrakqIkJCQgITJkwgMjKS7t2r7qwxPDyKnj3fYvv2z9m48T3S078lO3sTnTs/R0JCapXlML7l1zNWEXlCRNaJyFoRmSQiMSJSS0T+LSJbnN8XeS0/QkS2isgmEbnBa3pHEVnjzHtNTv13b0wlFRcX89hjj/GHP/zh9LRevXpVaVE9JSwsgubNB9Cr1zvEx6dy9Og2Zs/+NQcOLK3yLMY3/FZYRaQB8BjQSVXbAeHAAOBpYIaqtgBmOO8RkTbO/LZAP+AtEQl3NjcWGAq0cH76+Su3CX3FxcU8+OCDvP7667z++uvs2LHD7UgAJCU14+qr3yUlpRdFRcfYvn2K25FMJfm7jTUCiBWRCCAO2AvcBox35o8H+juvbwMmq+pJVd0BbAW6iEgKkKSqC9TT4j/Bax1jKqSoqIghQ4YwYcIE4uPjmT59Ok2aNHE71mmRkQl07vw32rV7jI4d/9vtOKaS/FZYVXUP8DKQDmQCOar6L6CuqmY6y2QCdZxVGgC7vTaR4Uxr4LwuPf0MIjJURJaKyNIDBw748nBMCFBVHn74YT7++GMSEhL47rvvuOaaa9yOdQaRMJo1u5OIiFgACgsLmTlzpsupTEX4syngIjxnoU2A+kC8iNx7rlXKmKbnmH7mRNVxqtpJVTvVrm1dB83PjRkzhvfff5/Y2FimT59Oz5493Y50XqrF3H777fTp04fPP//c7TimnPzZFNAH2KGqB1S1EJgCXAnsd77e4/zOcpbPABp6rZ+Kp+kgw3lderox5Zafn89bb70FwIQJE7jyyitdTlQ+IuF07NiRkpISBg4cyH/+8x+3I5ly8GdhTQe6iUiccxX/OmADMA24z1nmPuCfzutpwAARiRaRJnguUi12mgtyRaSbs50hXusYUy6xsbEsWLCA8ePH88tf/tLtOBUycuRInnjiCQoLC7nzzjvZvHmz25HMefizjXUR8DmwHFjj7Gsc8ALQV0S2AH2d96jqOuBTYD3wPTBcVYudzQ0D3sNzQWsb8J2/cpvQVadOHYYMGeJ2jAoTEV5++WVuu+02srOzufXWW8nOznY7ljkH8e5aF0o6deqkS5fafYDV3fPPP8/Ro0d59tlniY6O9ss++vf/0S/bPWXq1B4AHDt2jCuvvJI1a9bQr18/vvnmG8LCrPNkVRORZara6VzL2KdiQtaSJUsYOXIko0ePZsmSJW7HuWAJCQlMmzaN5ORkioqKyM3NdTuSOQvr0mpCUmFhIQ888ADFxcU88cQT9OjRw+1IPpGWlsaCBQto2rSpna0GMCusJiS99dZbrFu3jqZNm/L888+7Hcenmjdvfrr5QbWY4uKC0/e8+sqp5gdTOfZfngk5WVlZjBw5EoBXXnmF2FjfFp1AkZ+fxbx5v2PVqv9xO4opxQqrCTkjRowgJyeHfv36ccstt7gdx2+Ki09y5MgGMjKms3//ArfjGC9WWE1IKSoqYv/+/URGRvLqq68SygOhJSQ0pHXrhwBYteplCguPu5zInGKF1YSUiIgIvvrqK1avXk3Lli3djuN3TZveSc2arcjPz2L9+rfdjmMcVlhNyBERWrVq5XaMKhEWFsEVVzyNSAQ7d07l8OF1bkcyWGE1IWTkyJEsX77c7RhVLimpGc2bDwBgzZr/RbXE5UTGCqsJCXPmzOHZZ5+ld+/eHD9e/doaW7YcTHx8KvXq9bDCGgDsPlYTEk7dXvXEE08QHx/vcpqqFxERR+/eEwkLs3/SgcA+BRP0lixZwqxZs0hKSuLxxx+vsr77gca7qKoW89OTjUxVs6YAE/TGjBkDwNChQ6lRo4bLady3Y8eX/Oc/95CXt9/tKNWWFVYT1NLT0/nss8+IiIjgsccecztOQDh0aBV5eZls2vSh21GqLSusJqi99tprFBcXc9ddd9GwYcPzr1ANtGr1K0TCSU//jtzcnW7HqZasjdUEtUceeYTi4mIGDx7sdpSAkZDQkMaNb2bnzn+yYcN7dOnynNuRqh07YzVBLS0tjVdeeYUOHTq4HSWgtGx5P2FhUWRmzubIkY1ux6l2rLCaoKSqhOrTL3whNjaZpk1vB2Dz5gkup6l+rLCaoDR37lwuv/xyPvroI7ejBKxmzQYQFhbFsWPpFBXlux2nWrE2VhOUxo0bx+rVq9mwYYPbUQJWTMzF9Ow5lho1mtk9rVXMCqsJOocOHeLzzz9HRPj1r3/tdpyAVrNm6I/wFYisKcAEnfHjx3Py5EluuOEG0tLS3I4TFPLzs9i9+19ux6g27IzVBBVVZdy4cQD85je/cTlNcCgoyGXGjIGUlBSRnHwFsbG13Y4U8uyM1QSVhQsXsmnTJurVq8fNN9/sdpygEBWVSN26V6JazI4dX7gdp1qwwmqCyqRJkwC49957iYiwL1zl1azZ3QDs3PlPioryXE4T+uxvpgkqL7zwAt27d6djx45uRwkqtWq1pVatdhw+vJb09G9p2vSXbkcKaXbGaoJKXFwc99xzT7V4npWvNWvmecrAtm2foVrscprQZoXVBI2ioiK3IwS1lJQexMXVJy9vL5mZ/h2ztrqzwmqCwsGDB6lXrx5Dhw61rqyVJBJO8+YDSEm5mri4em7HCWnWxmqCwieffMKhQ4fYvXs3IuJ2nKDVpMl/0aTJf7kdI+RZYTV+54tHpcyd+w4AWVmdz9heoD4qxVRf1hRgAl5+/n4OH15DWFgU9epZEfWFrKxFLFny3xQWVr8n2lYFK6wm4O3ZMxOAunW7ExkZ53Ka0LB580T27p3F7t3T3Y4SkqywmoC3Z88PAKSmXudyktDRpIlnrNadO7+0i4F+YIXVBLS8vP1kZ28gPDyWOnW6ux0nZKSk9CI6uha5uTs5fHiN23FCjhVWE9BiY+tw7bV/54or/khERIzbcUJGWFgEjRrdCMCuXV+5nCb0WGE1AU1ESEpqRoMG1gzga40b3wJ4mloKCnJdThNarLCagGVtf/4VH1+f2rU7UVJSQEaGjdXqS34trCJSU0Q+F5GNIrJBRLqLSC0R+beIbHF+X+S1/AgR2Soim0TkBq/pHUVkjTPvNbE7xKuFrVs/ZvbsoezbN9/tKCGrWbO7adXqIVJSerkdJaT4+4z1VeB7VW0FXAZsAJ4GZqhqC2CG8x4RaQMMANoC/YC35KcH9YwFhgItnJ9+fs5tAsDevbPIzt5gA4b4Ud263bjkkvts8Gsf81thFZEkoBfwPoCqFqhqNnAbMN5ZbDzQ33l9GzBZVU+q6g5gK9BFRFKAJFVdoJ7vhhO81jEhKi9vH9nZG527Abq4HceYCvHnGWtT4ADwoYisEJH3RCQeqKuqmQDO7zrO8g2A3V7rZzjTGjivS08/g4gMFZGlIrL0wIEDvj0aU6UyM+cAnjOq8PBol9OENtUSNm2awKxZD9pjsn3En4U1AugAjFXVK4DjOF/7z6KsdlM9x/QzJ6qOU9VOqtqpdm37ahPM9u6dDUD9+le7nCT0iYSxf/98cnK2nP5zNxfGn4U1A8hQ1UXO+8/xFNr9ztd7nN9ZXss39Fo/FdjrTE8tY7oJUSdOHDo9NoB1Cqgap+5p3b37W5eThAa/FVZV3QfsFpFLnEnXAeuBacB9zrT7gH86r6cBA0QkWkSa4LlItdhpLsgVkW7O3QBDvNYxISgzcy6g1KnT2cYGqCINGlxHeHg0Bw+u4PhxO2+5UP4eNvBR4B8iEgVsBx7AU8w/FZFfAenAnQCquk5EPsVTfIuA4frT5eBhwN+BWOA758eEqNTUvkRGJhAdXcvtKNVGZGQ8KSnXkJExnfT0b4G73I4U1PxaWFV1JdCpjFlldqNR1eeB58uYvhRo59NwJmBFRsaTmtrH7RjVTqNGN5KRMZ3du7+juLiY8PDw869kymQ9r4wxACQnX05cXAr5+VnMmjXL7ThBzZ4gYALKqlX/Q0lJAS1bDiE+vsy76oyfiITRps0wwsOj6NXLemJdCCusJmAUFxewe/d0iovzueSS+92OUy01aHAtAJGRkS4nCW7WFGACxsGDKyguzicpqTlxcSlux6n2bBCcyrPCagLGvn1zAUhJ6elykuotL28/d999NzfddJPbUYKWNQWYgKBa4ty/ij0w0GWRkfFMmzaNEydOsHPnTtLS0tyOFHTsjNUEhCNHNnDy5GFiY+tSo0YLt+NUa5GRCfTv3x+AiRMnuhsmSFlhNQFh374fAahX7ypsuF333Xefp3PkxIkTra21EqwpwASE1NQ+iIRTt66NDRAI+vTpQ7169diyZQuLFi2iW7dubkcKKnbGagJCUlIzWrd+iFq12rodxQAREREMGjQIsOaAyrDCaowp05AhQwD45JNPKCoqcjlNcLHCaly3du3rbN06icJCe1JoILn00kt54403WLJkCRER1mpYEfanZVxVUJDL9u1fANCo0c0upzGlDR8+3O0IQcnOWI2r9u9fgGoxF198GVFRiW7HMWehqtYcUAFWWI2r9u2bB1ingED22Wef0bZtW8aNG+d2lKBhhdW4pqSkkKyshYDn/lUTmAoLC9mwYQMTJkxwO0rQsMJqXHPw4EqKivJITGxKfHx9t+OYs+jfvz8JCQksWrSITZs2uR0nKFhhNa7x7m1lAldcXBx33nknYPe0lpcVVuOa2rU7kpJyNSkpNqhyoDt1T+vEiRMpKSlxOU3gO2dhFZE7nd9NqiaOqU5SUnrRpctzXHRRK7ejmPPo1asXjRo1Ij09nTlz5rgdJ+Cd74x1hPP7C38HMcYErrCwMAYPHgzAp59+6nKawHe+DgKHRGQm0EREppWeqaq3+ieWCXWbN0+gRo2W1K7dibAw66cSDB566CE6duxoA2CXw/n+Rt8EdAAmAv/j/zimOsjPP8CGDe8SHh7NL37xDdYBMDikpaXZoNfldM6/0apaACwUkStV9UAVZTIh7lSngNq1OxMeHu1yGlMZeXl5xMXFuR0jYJ2zsIrIV4A6r8+Yb00BpjJ+us3KelsFG1Xl17/+NZMmTWL9+vU0btzY7UgB6XwXr17G0wSwA8gH3nV+jgFr/RvNhKLCwuMcOLAMCKNevSvdjmMqSETIy8sjLy/P7mk9h3MWVlWdraqzgStU9W5V/cr5GQjY6YapsKysRagWUatWO6KjL3I7jqmEU49tGT9+vD225SzK20Ggtog0PfXGua+1tn8imVB2qhkgJcX+Xw5Wffr0oX79+mzdupV58+a5HScglbewPgHMEpFZzu1XM4HH/ZbKhKyEhMbExze09tUgFh4efron1vvvv+9ymsBU3sI6C3gHOILnYtY7wGw/ZTIh7JJL7qNPn49JSGjodhRzAR588EHA01ng6NGjLqcJPOUtrBOAJsBrwN+c19ZybUw11aJFC3r16kVeXh4zZ850O07AKe+d2Zeo6mVe72eKyCp/BDKhSVXZtesr6tTpSlxcXbfjGB945ZVXSExMpEWLFm5HCTjlLawrRKSbqi4EEJGugLVam3LLydnMqlUvEROTzPXXTynzvmgTXDp06OB2hIBV3qaArsB8EdkpIjuBBcDVIrJGRFb7LZ0JGZmZnhGR6tXraUU1xKgqe/bscTtGQCnvGWs/v6YwIW/vXs+1zvr1bezVUHLo0CGuvfZa9u7dS0ZGBjExMW5HCgjlOmNV1V3n+vF3SBPccnN3cuzYLiIjE7n44svdjmN8qFatWkRERHDo0CGmTJnidpyAYU8QMH73UzPAVTZEYIgRER5++GEA3n77bZfTBA6/F1YRCReRFSLytfO+loj8W0S2OL8v8lp2hIhsFZFNInKD1/SOTnvuVhF5TayRLqhkZs4FICXlapeTGH+45557SExMZO7cuaxbt87tOAGhKs5Yfwds8Hr/NDBDVVsAM5z3iEgbYADQFk+b7lsiEu6sMxYYCrRwfqzNN0gUFBQQHV2LyMgk6tTp7HYc4weJiYnce++9ALzzzjsupwkMfi2sIpKKZ7Ds97wm3waMd16PB/p7TZ+sqidVdQewFegiIilAkqouUM+IDxO81jEBLioqim7dXqRfv6k29moI+81vfgPAhAkTOH78uMtp3OfvM9b/BZ4CvB/rWFdVMwGc33Wc6Q2A3V7LZTjTGjivS083QSQsLNLtCMaPLrvsMrp3705ERAQbNmw4/wohzm9XEkTkZiBLVZeJyDXlWaWMaXqO6WXtcyieJgMaNWpUvqDGb7Kysli9ejUlJRF20aoamDBhAg0aNCA2NtbtKK7z5xnrVcCtToeCyUBvEfkI2O98vcf5neUsnwF4j8yRCux1pqeWMf0MqjpOVTupaqfatW1UQ7dNnjyZvn37smLFC25HMVWgefPmVlQdfiusqjpCVVNVNQ3PRakfVPVeYBpwn7PYfcA/ndfTgAEiEu2M99oCWOw0F+SKSDfnboAhXuuYAHbqMcl163ZzOYmpSrm5uXz99ddux3CVG9/PXgA+FZFfAenAnQCquk5EPgXWA0XAcFUtdtYZBvwdiAW+c35MAMvIyGDevHnExMTYI1iqkfz8fJo2bcrhw4fZtm1btX2qa5V0EFDVWap6s/P6kKpep6otnN+HvZZ7XlWbqeolqvqd1/SlqtrOmfeI2vMgAt4XX3wBwI033khEhD3Ns7qIjY2lX79+lJSU8Oabb7odxzXW88r4xalmgLvuusvlJKaqPfbYYwC899575ObmupzGHVZYjc+lp6czf/58YmJiuOmmm9yOY6pY586d6dGjB9nZ2bz77rtux3GFFVbjc7t27aJJkybceuutJCQkuB3HuOCpp54CYMyYMRQUFLicpupZYTU+17NnT7Zt28a4cePcjmJcctNNN9GmTRv27NnDxx9/7HacKmeF1fiFiFCjRg23YxiXhIWF8dRTT9G1a1caNqx+D460wmp86scff2Tv3jL7b5hqZvDgwSxYsIDrrrvO7ShVzgqr8ZmSkhIGDRpEamoqq1bZsyaru7CwsGr7GB4rrMZn5s2bR3p6OqmpqbRv397tOCZAbNmyhSFDhvDNN9+4HaXK2MgYxmc++ugjAAYOHEhYmP2fbTy+/vprJk6cyPr167nxxhurxVms/e03PpGXl8fkyZMBT9uaMac8/PDD1KtXj2XLlvHVV1+5HadKWGE1PjFlyhSOHj1K165dadu2rdtxTACJjY3l6aefBmDkyJGUlJScZ43gZ4XV+MQHH3wAwIMPPuhyEhOIhg4dSv369Vm5cuXp7s6hzAqruWCFhYXUqFGDpKQk7r77brfjmAAUGxvLqFGjABgxYgQnT550N5CfWWE1FywyMpIvv/ySPXv2WKcAc1YPPPAAbdq0IT09nTlz5rgdx6/srgDjMzYugDmXiIgIPvjgAxITE2nTpo3bcfzKCqu5IIsXL+bIkSP07dvXbrEy59W1a1e3I1QJ+5dgLsioUaPo168fb7/9tttRTBBRVaZMmcL69evdjuIXVlhNpW3bto3vv/+e6OhoG9DaVMgbb7zBHXfcwfDhwwnFB4JYYTWVNnbsWFSVAQMGkJyc7HYcE0QGDRpEcnIys2bNYtKkSW7H8TkrrKZS8vLyTt+7+sgjj7icxgSbWrVq8eKLLwLw5JNPkpOT43Ii37LCaipl8uTJHDlyhC5dutCpUye345ggdP/999O9e3f27dvHH//4R7fj+JQVVlNhqsqrr74KwPDhw11OY4JVWFgY48aNIzIyknfeeYcffvjB7Ug+Y7dbmQorLCzk7rvvJiwszHpaVRP9+//ot203azaEjRvf58UXX6R3795+3x/A1Kk9/Lp9O2M1FRYVFcWf/vQnli9fTnR0tNtxTJBr0eJeXnrpJaZOnep2FJ+xM1ZTadVhXE3jf2FhEfz+9793O4ZP2RmrqZDf/va3jB49mtzcXLejmBB09OhRfve733HyZLbbUS6InbGactu4cSNvv/02kZGRDB48mMTERLcjmRAzbNgwPv74Y1JSVtK583NB+63IzlhNuY0aNQpV5f777yclJcXtOCYEPf/88yQlJZGZOYddu6a5HafSrLCaclm5ciWffPIJUVFR/PnPf3Y7jglRaWlpjB07FoA1a14jJ2eLy4kqxwqrKZdTxfS3v/0tDRs2dDmNCWUDBw6kceNbKCkpYMmSv1BYeNztSBVmhdWc1/z58/nmm2+Ij49nxIgRbscx1UD79r8jKak5x49nsHLl6KAbqMUKqzmvzz77DIAnnniCOnXquJzGVAfh4dF07vxXwsNjiYpKRLXY7UgVYncFmPMaM2YMnTt35sYbb3Q7iqlGEhIa0bv3ROLi6rodpcKssJrzEhEGDhzodgxTDXkX1cLCXIqKThIbG/hDVFpTgDmriRMnsnXrVrdjGMOxY+nMnj2UxYtHUFwc+E94tcJqyrRx40Z+9atf0b59ezIzM92OY6q5qKgaqBaTnb0xKC5mWWE1ZygqKuKBBx6gsLCQQYMGWWcA47qoqBp06fJ/CQ+PISPjX2zb9onbkc7JCqs5w+jRo1m4cCENGjTg5ZdfdjuOMQDUqNGcDh2eAWDdurFkZS1yOdHZ+a2wikhDEZkpIhtEZJ2I/M6ZXktE/i0iW5zfF3mtM0JEtorIJhG5wWt6RxFZ48x7TYK1A3EQWLlyJaNGjQLgww8/pGbNmq7mMcZb/frX0LLlfUAJS5aMIjc33e1IZfLnGWsR8KSqtga6AcNFpA3wNDBDVVsAM5z3OPMGAG2BfsBbIhLubGssMBRo4fz082Puauv48ePce++9FBYWMnz4cPr27et2JGPO0KrVg6Sk9KKo6Bj7989zO06Z/Ha7lapmApnO61wR2QA0AG4DrnEWGw/MAv7oTJ+sqieBHSKyFegiIjuBJFVdACAiE4D+wHf+yl5dLV++nO3bt9OyZcvTD3ozJtCIhNGhwzNkZfWlfv1r3I5Tpiq5j1VE0oArgEVAXafooqqZInKqK08DYKHXahnOtELndenpZe1nKJ4zWxo1auTDI6geevbsybx584iJiSE+Pt7tOMacVURE3M+KalHRCSIiYtwLVIrfL16JSALwBfC4qh4916JlTNNzTD9zouo4Ve2kqp1q165d8bDV1JEjR06/vuKKK2jdurWLaYypmJycrcycOYT09G/djnKaXwuriETiKar/UNUpzuT9IpLizE8BspzpGYD3sEmpwF5nemoZ040PzJkzh7S0tNPjARgTbI4cWU9eXiarVr3M4cNr3Y4D+PeuAAHeBzao6hivWdOA+5zX9wH/9Jo+QESiRaQJnotUi51mg1wR6eZsc4jXOuYCLFmyhJtvvpmjR48yY8YMt+MYUylpabfSpMntlJQUsnjxM+TnZ51/JT/z5xnrVcBgoLeIrHR+bgReAPqKyBagr/MeVV0HfAqsB74HhutPQ9oMA94DtgLbsAtXF2z+/PnccMMN5ObmMmDAAN588023IxlTae3aPUpycgdOnjzMokV/cr3bqz/vCviRsttHAa47yzrPA8+XMX0p0M536aq3L7/8koEDB3LixAluu+02JkyYQHh4+PlXNCZAhYVF0KnTX5kzZyg5OZtYufJFOnT4b9eemWU9r6qZ8ePHc8cdd3DixAmGDh3K559/TmRkpNuxjLlg0dE16dr1/xEeHktGxr85fHi1a1ls2MBqpkOHDkRHR/PMM8/wzDPPBO1TMI0pS1JSMzp0+DMlJQVcfPFlruWwwloNzJkzh549eyIitG/fns2bN9tzq0zIql+/l9sRrLAGgv79f/TLdo8f38vata+xb988/vGPf5werPrRR3cBu/yyT4CpU3v4bdvGVMSRIxvZuPFdOnX6G5GRcVW2XyusISg//wCbN09k166vUC0iKSmJgoICt2MZU6VUlVWrXiYnZxPLlz9Hly7PIVI1l5WssIaQ3NydbNv2Cbt3/4uSkgJASE29gcWLP7QxVU21IyJ06jSSOXOGsm/fXDZtGk+rVg9Uyb7troAQcujQKnbt+pqSkgLq17+Ga68dT8eOf7aiaqqthISGdOw4Cghj06YPyMycUyX7tTPWIFRcXMCBA0vYs2cGMTF1aNv2YQBSU2/g2LHdpKXdSkKCDUJjDEDdul1p02Yo69e/zbJlz9Gr1zuAf68DWGENEoWFx9i/fyH79s1l//6FFBXlARAZmUTr1g8RFhZBREQM7do94nJSYwJP8+YDycnZyp49/2Hx4j9x8uRdREdH+21/VliDwK5dX7Nq1cv81MMXkpKa06BBb+rXv5awMPsYjTkXEeHyy//IiRMHaNr0Tr8WVbDCGlBKSoo4fHgt+/fPIzGxGY0aeR6UkJjYBICLL76MlJSe1KvXk/j4+m5GNSboRETEcNVVr1dJpxgrrC7Ly8sjM3MumZlz2LdvPoWFniFrL7748tOF9aKLWtOv3zSiopLcjGpM0PMuqgsWLGDXrl0MGDDA5/uxuwJcNHr0aGrXrs3ixX9i9+7vKSw8Snx8Ks2aDaB164dOLycSZkXVGB/asmULt99+O8nJyX7Zvp2xluKvXlAlJYXs37+I99+/mXbtPAN11alTh7y8PGrWbE1KSk9SUnqSmJjml/0bY37SokULVq9ejb+eNGKF1Y9UlezsjaSnf8fevT9QUJDDm2/uYuzYsQDccccd9O7dm8ceC8xH+BoTyvz5+CYrrH5w8mQ2u3d/T3r6t+Tm7jg9PTGxCZdeeqnX+0QSExMBK6zGhBIrrH6wefN4tm//HICoqJo0bHg9DRv2IympOcOG9XQ5nTHG36ywXqCiohNkZPyLmJiLqVfvKgAaN76FY8cyaNz4FurV605YmA0kbUx1YoW1kvLzs9i+fQq7dk2jsDCXmjVbny6sSUlN6d79JZcTGmPcYoW1gnJytrJ16yT27JlxuifURRe1oUmTO1BVG5HfGGOFtSL27p3NkiV/dt6FUb9+b5o1u4tatdq6mssYE1issJ6DqpKXl3m6+2idOl2Ija1HSkoPmjW7i7g4G47PGHMmK6xlUNXTA+Pm5e3j+us/IyIijoiIWPr0mWSDnhhjzskqhJeSkhL27p3Npk1/5+jRrQBER9ciN3cXF13UGsCKqjHmvKxK4DlDnTp1KqNGjWL1as+zyGNikmnefBBpabcQHu7fIcaMMaHFCqvjueeeY/Xq1cTEJNOixWAaN77JCqoxplKssOIZSuyFF15g48aN/OtfbaygGmMuiA0b6Ojbty+PPvqoFVVjzAWzwmqMMT5mhdUYY3zMCqsxxviYFVZjjPExK6zGGONjVliNMcbHrLAaY4yPWWE1xhgfs8JqjDE+ZoXVGGN8LGgKq4j0E5FNIrJVRJ52O48xxpxNUBRWEQkH3gR+AbQB7hGRNu6mMsaYsgVFYQW6AFtVdbuqFgCTgdtczmSMMWUKlmEDGwC7vd5nAF1LLyQiQ4GhzttjIrKpCrKdSzJw0HuCGw9x9fM+XT/GKtjfz44xBI8PqtkxXuD+Gp9vgWAprGX9MegZE1THAeP8H6d8RGSpqnZyO4c/2TGGBjtG3wqWpoAMoKHX+1Rgr0tZjDHmnIKlsC4BWohIExGJAgYA01zOZIwxZQqKpgBVLRKRR4DpQDjwgaquczlWeQRMs4Qf2TGGBjtGHxLVM5oqjTHGXIBgaQowxpigYYXVGGN8zAprKefrOiserznzV4tIB695H4hIloisLbXOnSKyTkRKRKST1/Q0EckXkZXOz9te8zqKyBpnP6+J+O5OvwA6xllOjlPz6gTjMTrzLhWRBc78NSIS40wPic/xPMcYEp+jiAzyOoaVzvzLnXkV+xxV1X6cHzwXxrYBTYEoYBXQptQyNwLf4bm3thuwyGteL6ADsLbUOq2BS4BZQCev6Wmll/Watxjo7uznO+AXIXiMP1s2iD/HCGA1cJnz/mIgPMQ+x3MdY0h8jqWWaQ9s93pfoc/Rzlh/rjxdZ28DJqjHQqCmiKQAqOoc4HDpjarqBlUtdy8wZ3tJqrpAPZ/qBKB/pY7oTAFxjH5W1cd4PbBaVVc5yx1S1eIQ+xzLPEYfHcvZuPl39R5gElTu36MV1p8rq+tsg0osUxFNRGSFiMwWkZ5e+8jw4T68BcoxnvKh87Xrv334Nbmqj7EloCIyXUSWi8hTXvsIlc/xbMd4Sih8jt7uximsVOJzDIr7WKtQebrOlqt7bTllAo1U9ZCIdASmikhbH++jtIA4RlU9CgxS1T0ikgh8AQzGczZwoar6GCOAHkBnIA+YISLLgKM+3EdpAXGMqjqD0PkcPRsU6QrkqeqpttkK78POWH+uPF1nfda9VlVPquoh5/UyPO1JLZ19pPpiH2UIlGNEVfc4v3OBj/F89fOFKj1GZ1uzVfWgquYB3+Jp2wuZz5GzH2MofY6nDOCns9VT+6jQ52iF9efK03V2GjDEuRrZDchR1czK7ExEaotnrFlEpCnQAk+DeSaQKyLdnK9VQ4B/VvKYSguIYxSRCBFJdqZHAjcDa8++pQqp0mPE0yPwUhGJE5EI4GpgfSh9jpzlGEPsc0REwoA78bTnAlCpz7GiV+pC/QfPVcbNeM6snnGmPQw87LwWPINubwPW8PMrp5PwfPUtxPO/3K+c6f/lvD8J7AemO9PvANbhudq5HLjFa1ud8PwF3Qa8gdNLLlSOEYgHluG50rwOeBXnKnOwHaMz717nONYCo0PtczzbMYbg53gNsLCMHBX6HK1LqzHG+Jg1BRhjjI9ZYTXGGB+zwmqMMT5mhdUYY3zMCqsxxviYFVZTbiJy7ALWvVWc0YlEpL+ItPGa96yI9PFFxlL7nCWlRmiqzDLl3NfjIhLn9f5bEal5odsttY9rRCRHRL49z3IzReSYL47LVI4VVlMlVHWaqr7gvO0PtPGa9xdV/Y8rwXznceB0YVXVG1U12w/7mauqN55rAVW9Fljqh32bcrLCaipMRMJE5C3xjGn5tXN29ktn3k4R+aszUMcaEWnlTL9fRN4QkSuBW4GXnEE7monI373Wv048A7asEc94mtHn2m6pXLEiMlk843J+AsR6zbtePGOJLheRz0QkoYz1x4rIUue4/uqV50uvZfqKyJRS6z0G1AdmishMr7zJ4hmPdqOIvCcia0XkHyLSR0TmicgWEeniLB/vHO8S5/hLj+JU1ueQIiJznD/HtXLmADfGJVZYTWXcjmec1fbAQ3jGqfR2UFU7AGOB33vPUNX5eLoh/kFVL1fVbafmiWfg5L8Dd6tqezwDfwwrz3Ydw/AMnnEp8DzQ0dluMvBnoI+z/lLg/5Sx/jPqee78pcDVInIp8APQWkRqO8s8AHxY6phew9N3/FrnbLG05nh6JF0KtAIG4hnQ5PfAn07tG/hBVTsD1+L5jye+jG15G4in19DlwGXAyvMsb6qIFVZTGT2Az1S1RFX3ATNLzT91RrcMTwEur0uAHaq62Xk/Hs9gxeXdbi/gIwBVXY2nmyV4BkBuA8wTkZXAfUDjMta/S0SWAyuAtngGVVZgInCv02baHc9AxxWxQ1XXqGoJnm6fM5ztrvE6juuBp518s4AYoNF5trsEeEBERgHt1TMIigkANmygqYzzjbd50vldTMX+jvliu2X10Rbg36p6z1l3LNIEzxlkZ1U9IiJ/x1PcwHOG+hVwAs9/KEXnyXm23AAlXu9L+Ok4BLhDKzBYuKrOEZFewE3ARBF5SVV9MVyfuUB2xmoq40fgDqettS6egSsqIhdILGP6RiBNRJo77wcDsyuw3TnAIAARaYfnqzfAQuCqU9sVzwhNLUutmwQcB3KcY/rFqRmquhfPV/0/42mqqMgxldd04FFn9CRE5IrzrSAijYEsVX0XeB9nGD/jPiuspjK+wDM60FrgHWARkFOB9ScDf3Au0jQ7NVFVT+Bpw/xMRNbgOaN7+yzbKMtYIEFEVgNP4XlOEap6ALgfmOTMW4inrfM09TxyZAWer+ofAPNKbfsfwG5VXX+WfY8Dvjt18aoS/gZEAqvF8/C7v5VjnWuAlSKyAs8oYq9Wct/Gx2x0K1MpIpKgqsdE5GI8Bewqp701JInIG8AKVX3fxQzXAL9X1ZvLsewsZ1m77coF1sZqKutr52JOFPC3EC+qy/A0EzzpcpQCoJ2IfHuue1mds+ameMYhNS6wM1ZjjPExa2M1xhgfs8JqjDE+ZoXVGGN8zAqrMcb4mBVWY4zxsf8PcjIhOYko5GQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=[5,4])\n",
    "n, bins, patches = plt.hist(x=IgnTime, bins='auto', color='#0504aa',\n",
    "                            alpha=0.7, rwidth=0.85,density=True)\n",
    "plt.plot(x, gkde.evaluate(x), linestyle='dashed', c='black', lw=2)\n",
    "plt.xlabel('Ignition delay time [s]')\n",
    "plt.ylabel('pdf')\n",
    "# plt.ylim([0,4])\n",
    "plt.savefig('GRi3IgniDelayPdf.pdf', bbox_inches = 'tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "del(tchem)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "pytchem.finalize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
